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PREFACE 


This  Memorandum  is  one  in  a  continuing  series  of  RAND 
publications  dealing  with  theoretical  computational  ques¬ 
tions  arising  from  the  RAND  program  of  research  in  biology 
and  physiology.  The  Memorandum  contributes  to  our  ability 
to  apply  computer  technology  to  the  analysis  of  complex 
chemical  systems  by  considering  the  "chemical  equilibrium 
problem,"  the  problem  of  determining  the  distribution  of 
chemical  species  that  minimizes  the  free  energy  of  a  system 
while  conserving  the  mass  of  each  of  the  chemical  elements. 

Solutions  to  the  chemical  equilibrium  problem  pub¬ 
lished  up  to  this  t  .me  [4,5]  apply  to  those  problems  for 
which  an  estimate  of  the  solution  exists.  This  Memorandum 
considers  a  problem  for  which  no  estimated,  solution  exists 
and  solves  that  problem  with  the  maximum  precision  now 
available. 

The  mathematical  aspects  of  this  Memorandum  should 
also  be  of  interest  in  other  fields  where  computational 
analyses  of  complex  chemical  systems  are  under  considera¬ 
tion,  e.g.,  in  studies  of  rocket  propulsion  systems, 
planetary  atmospheres,  re-entry  problems,  etc. 


SUMMARY 


In  physical  chemistry,  the  "chemical  equilibrium 
problem"  is  the  problem  of  determining  the  distribution 
of  chemical  species  that  minimizes  the  free  energy  of  a 
system  while  conserving  the  mass  of  each  of  the  chemical 
elements.  The  reactions  occurring  within  the  chemical 
system  may  be  quite  complex.  However,  in  a  great  number 
of  cases,  the  mathematical  statement  of  the  problem  can 
be  simplified  to  a  particular  mathematical  form  [7,8] 
involving  the  minimization  of  a  nonlinear  objective  func¬ 
tion  over  a  set  of  linear  constraints. 

This  Memorandum  presents  the  numerical  solution  of 
the  chemical  equilibrium  problem  by  describing  methods 
for  starting  the  solution  when  an  initial  estimate  is  not 
available,  and  for  improving  an  initial  estimate  to  make 
it  feasible.  It  presents  a  first-order  method  and  a 
second-order  method  for  solving  the  chemical  equilibrium 
problem  in  the  context  of  the  linear-logarithmic  program¬ 
ming  problem  [4]  and  provides  convergence  criteria  for 
the  majority  of  problems  of  this  type  that  are  likely  to 
be  attempted. 
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FOREWORD 


In  deciding  between  the  languages  of  mathematics  and 

physical  chemistry,  we  have  chosen  in  this  Memorandum  to 

use  that  of  mathematics.  The  disadvantage  of  this  choice 

is  that  the  physical  chemist  may  experience  some  difficulty 

in  immediately  identifying  certain  concepts.  The  advantage 

is  that  mathematical  language  divorces  the  methods  from 

the  physical  assumptions  involved  in  constructing  a  mathe- 

* 

maticai.  model  of  a  physical  system.  The  mathematical 
methods  are,  hence,  free  to  transcend  their  specific 
chemical  applications. 

The  methods  given  here  do  not  solve  every  problem  that 
is  specified  in  the  given  mathematical  form.  The  solution  of 
a  problem  in  which  some  phase  vanishes  (a  degenerate  problem) 
requires  further  work.  Some  work  has  been  done  on  particular 
degenerate  systems  [13],  but  the  accurate  numerical  solution 
of  a  large  general  system  of  this  type  has  yet  to  be  accom¬ 
plished.  Until  recently,  a  skilled  physical  chemist  could 
intuitively  eliminate  the  degeneracies  of  his  model  and 

* 

The  reader  is  referred  to  other  works  for  the  pro¬ 
cedure  of  constructing  the  mathematical  models  of  bio¬ 
chemical  systems  [9-12]. 
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obviate  the  need  for  solving  a  degenerate  system.  But, 
as  problems  grow,  eliminating  degeneracy  becomes  increasingly 
difficult.  Frequently,  the  point  at  which  the  problem  be¬ 
comes  too  large  for  the  physical  chemist  to  decide  whether 
or  not  to  Include  a  phase  coincides  with  the  point  at  which 
the  problem  becomes  numerically  unwieldy.  Hopefully,  the 
future  will  eliminate  these  difficulties. 

Statements  about,  convergence  and  convergence  tests 
exist,  unless  otherwise  indicated,  in  the  context  of  finite- 
accuracy  numerics.  Statements  of  this  kind  do  not  mean, 
in  the  absence  of  qualification,  that  no  problem  exists 
nor  that  no  machine  would  serve  as  a  counter  example. 

Rather  they  are  simply  descriptions  of  what  was  found  to 
occur  in  actual  practice. 

No  attempt  has  been  made  to  describe  those  methods 
which  were  tried  and  found  wanting.  The  methods  presented 
are  those  which  are  best  for  the  largest  number  of  cases. 

Finally,  it  should  be  pointed  out  that  although 
computing  time  was  a  factor,  it  was  consieered  secondary 
to  accuracy  of  results. 
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1 .  INTRODUCTION 


For  the  purposes  of  this  Memorandum,  the  chemical 
equilibrium  problem  is  merely  a  name  we  use  for  a  par¬ 
ticular  mathematical  programming  problcnn,  i.e.,  the  prob¬ 
lem  of  minimizing  a  particular  nonlinear  function  F(Xj^,X2, 
...,x^),  defined  below,  while  satisfying  the  linear  re¬ 
straints  or  constraints 

n 

y  a_  Xj  »  bi  i-l,2,3,...,m  (1.1) 

j  =  -'- 


with  x^  ^  0  for  j=l,2,...,n  and  given  constants. 

Assuming  that  the  equations  of  (1.1)  are  linearly  inde¬ 
pendent,  then  in  order  to  have  a  non- trivial  problem  it  can 
be  assumed  that  m<n.  The  variables  x, ,x„,...,x  can  be 
considered  components  of  a  vector  (x^^ ,X2 , .  . .  ,x^)  .  Solving 
the  chemical  equilibrium  problem  then  is  the  problem  of 
determining  this  vector.  The  variable  x^  will  be  referred 
to  as  the  "j  component";  also  the  numerical  value  of  x^ 
may  be  referred  to  as  the  "component"  rather  than  using 
the  perhaps  linguistically  correct  but  cumbersome  term 
"component  value." 
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The  components  are  partitioned  into  p  non-empty 
subsets  called  compartments .  Let  us  denote  these  compart¬ 
ments  by  <1 ),  <2 <p) .  Then  if  the  component  is 
t  h 

in  the  k  compartment,  we  will  say  j€<k),  where  each 
component  is  in  exactly  one  compartment.  The  number  of 
the  compartment  that  the  component  is  in  is  denoted 
by  [j].  Hence  j€ (k)  implies  [j]  =  k,  and  conversely. 

Each  compartment  has  associated  with  it  a  sum  defined  by 


I  • 

j€<k) 


(1.2) 


A  A  ^ ' 

The  component  fraction  x,  is  defined  by  x.  ■ 

J  J  ; 


whenever 


®[j]  ^  °- 


[j] 


The  objective  function  to  be  minimized  over  (1.1) 


is 


n 


F(Xj^,X2,  .  .  .  ,x^)  -  ^x^(Cj  +  log  -w^)  (1.3) 

j-i 


where  >  ^2  ’  '  *  ' ’^n  given  constants,  called  objective 
constants . 

When  an  X.  is  zero,  log  x.  is  undefined:  but  we  de- 
J  J 

fine  0  log  0  to  equal  0  so  that  we  may  evaluate  F  when 
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some  components  are  zero.  A  feasible  solution  to  the 
chemical  equilibrium  problem  is  defined  to  be  any  set  of 
non-negative  components  that  satisfies  (1.1).  The  problem 
is  said  to  be  feasible  if  it  has  feasible  solutions.  If 
no  feasible  solution  is  arbitrarily  large  in  any  component, 
the  feasible  problem  is  said  to  be  bounded  feasible;  all 
practical  problems  with  which  one  might  have  occasion  to 
deal  are  bounded  feasible. 

A  solution  or  optimal  solution  to  a  bounded  feasible 
problem  is  any  feasible  solution  in  which  F(Xj^ , .  .  .  ,x^) 
attains  the  minimum  value  over  all  feasible  solutions.  A 
problem  which  has  optimal  solutions  in  which  some  component 
is  zero  is  called  degenerate ,  and  a  bounded  feasible  prob¬ 
lem  in  which  the  components  in  any  optimal  solution  are 
all  strictly  positive  is  called  a  non-degenerate  problem. 

It  has  been  shown  [1,  Theorem  12.1]  that  a  non-degenerate 
problem  has  exactly  one  optimal  solution.  Hence,  we  may 
speak  of  the  solution  to  the  problem.  Furthermore,  it  has 

■k 

also  been  shown  for  the  non-degenerate  problem  that  the 
minimization  of  F  is  equivalent  to  the  existence  of  numbers 
TT  ,  7^2  >  •  •  •  » ^ni’  Lagrange  multipliers,  which  satisfy: 

Ref.  1,  p.  18. 
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m 

Yj  “  Cj  +  log  .  j-1,2,3, . .  .  ,n  (1.4) 

i=l 

In  the  following  sections  we  derive  conditions, 

analogous  to  (1.4),  which  are  useful  in  solving  the  problem. 

In  Sec.  2  we  are  interested  in  finding  a  solution  to  (1.1) 

with  all  X.  >  0.  A  set  of  x.  which  satisfies  these  con- 
J  J 

ditions  is  called  a  positive  feasible  solution.  If  (1.1) 

is  satisfied  with  x,  ^  0,  we  have  called  such  a  result  a 

J 

feasible  solution.  The  theory  of  linear  programming  gives 
us  methods  of  finding  feasible  solutions  to  problems  with 
linear  restraints.  In  Sec.  2,  we  use  a  linear  programming 
technique  to  find  a  positive  feasible  solution.  In  Sec.  4 
we  show  how  to  modify  the  initial  ositive  feasible  solu¬ 
tion  to  get  the  solution  to  the  problem. 
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2.  THE  INITIAL  SOLUTION 


The  algorithms  presented  in  the  following  sections 
require  an  initial  positive  feasible  solution  in  order  that 
the  procedure  for  solving  the  problem  can  be  initiated. 
Frequently,  an  individual  with  a  problem  to  solve  will  be 
able  to  give  a  rather  accurate  estimate  of  its  optimal 
solution.  This  estimate  may  be  the  exact  solution  of 
another  problem  which  differs  from  the  one  being  considered 
in  relatively  minor  ways. 


THE  PROJECTION  METHOD 

Let  us  suppose  that  such  is  the  case,  and  let  us  de¬ 
note  the  estimate  of  the  components  by  y^^ ,  y2  >  •  •  •  » •  These 
values,  substituting  y^  for  Xj  in  Eq.  (1.1),  will  not 
generally  satisfy  (1.1),  being  somewhat  in  error.  Let  us 
denote  these  errors  by  g^,g2i  •  •  .  that  is,  let 


■1  -  ‘’i  -  I 

j=l 


■  •  •  I  in 


(2.1) 


Then,  we  wish  to  find  corrections  to  y^  such  that,  denoting 
the  corrections  by  0 , ,  we  have 


^  ®  i-l,2,...,m 

j-1 
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or 


n 

=  I  *  i*l,2,...,ni  (2.2) 

j=l 

The  0.  must  also  be  chosen  such  that  y.  +0.  >  0,  for  all 
J  J  J 

j.  We  cannot  guarantee  this  condition,  but  we  can  attempt 

to  choose  small  values  for  0..  One  way  to  do  this  is  to 

J 

minimize 


I 


subject  to  (2.2),  where  w,  is  the  "weight"  or  relative 

J 

importance  of  minimizing  0^.  This  reduces  to  the  problem 
of  finding  Lagrange  multipliers  such  that 

with 

n  m  /  n 

I“j*j  -  I 'if 

j-1  i-1  \j=l 


we  have 


0  . 


j-l,2,...,n  (2.4) 
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Equation  (2.4)  becomes 


W.0. 
J  J 


m 


i-1 


(2.5) 


and  substituting  (2.5)  into  (2.2)  we  have 


g 


1 


w . 
J 


i“l , 2 , . . . ,m 


(2.6) 


The  terms 


I 


a 


j  =  l 


w. 

J 


can  be  immediately  evaluated;  let  us  denote'  these  terms 
by 


:i 


(2.7) 


Note  that  q  .  =  q.  .  Then,  (2.6)  becomes 

m 

gi  *  •  i=l,2,...,m  (2.8) 
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Equation  (2.8)  Is  a  set  of  m  simultaneous  equations  In 

the  m  unknowns,  ,  1T2  >  •  •  •  > These  equations  may  be 

solved  for  tt,  , tj_  , .  .  .  ,  n  ,  and  then  these  values  may  be  sub- 

stituted  in  (2.5)  to  get  9^^ ,02 » •  •  •  remains  the 

question  of  choosing  values  for  the  weighting  factors 

w..  In  tests  of  this  method,  it  has  been  found  that 
J 

using 


yields  satisfactory  results.  The  choice  of  the  weighting 
factors  depends,  to  some  extent,  on  the  available  com¬ 
puters.  Using  these  weighting  factors,  we  can  summarize 
the  computation  of  0.  in  the  following  three  equations: 


q 


>A 


n 


I 

j=l 


a. .a. .y . 
A  iJ  J 


m 


I 

/-I 


I  "ij'j 


m 


i-1 


i“l , 2 , .  .  .  ,m 
''^“1 , 2 ,  .  .  .  ,m 


(2.9) 


i-l,2,...,m  (2.10) 


j-1,2, 


,n 


(2.11) 


where 
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Xj  -  yj  +  9j  .  J-1,2 . n  (2.12) 

The  Xj  from  (2.12)  will  satisfy  (1.1).  However,  the 
Xj  need  not  all  be  strictly  positive.  If  any  x^  is  zero 
or  negative,  this  method  of  obtaining  the  initial  solution, 
which  we  shall  call  the  projection  method,  has  failed.  If 
the  projection  method  fails,  or  if  no  initial  estimate  is 
provided,  then  a  linear  programming  method  may  be  used. 

THE  LINEAR  PROGRAMMING  METHOD 

The  terminology  used  in  linear  programming  is  similar 
to  the  terminology  used  above  in  describing  the  chemical 
equilibrium  problem.  The  statement  of  a  linear  program¬ 
ming  problem  includes  a  set  of  linear  restraints 
n 

a^jXj  -  b^  i»l,2,...,m  (2.13) 

J-1 

together  with  a  set  of  constants  C^,C2,C^, . . . ,C^,  called 
costs .  A  feasible  solution  to  a  linear  programming  problem 
is  any  set  of  non-negative  Xj  such  that  (2.13)  is  satisfied. 
The  costs  are  used  to  form  the  following  expression,  L, 
which  is  called  the  objective  function 
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n 

L  =  ^  .  (2.14) 

j=l 

For  every  set  of  feasible  x,,  we  can  evaluate  L.  The  set 

J 

of  feasible  x.  for  which  L  has  the  minimum  value  that  it 
J 

can  have  with  any  set  of  feasible  x^ ,  is  called  a  solution 

of  the  linear  programming  problem.  A  problem  which  has 

sets  of  feasible  x^  is  called  a  feasible  problem,  and  a 

problem  in  which  there  are  no  sets  of  feasible  x.  is  called 

J 

an  infeasible  problem.  An  infeasible  problem  has  no  solu¬ 
tions,  while  a  feasible  problem  has  at  least  one  solution. 

In  this  discussion,  we  will  not  be  concerned  as  to  whether 
a  problem  has  moij  than  one  solution:  we  will  only  be 
concerned  with  finding  a  solution  to  the  problem.  Since 
the  means  of  finding  a  solution  to  a  linear  programming 
problem  has  been  the  subject  of  many  papers  and  books,  we 
will  not  give  an  actual  method  of  solving  the  linear  pro¬ 
gramming  problem  here.  The  reader  may  refer  to  Dantzig 
[2]  for  a  complete  discussion  of  the  problem. 

The  problem  of  finding  a  feasible  solution  to  a 
linear  programming  problem  is  itself  a  linear  programming 
problem-- that  is,  it  involves  finding  a  solution  to  the 
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problem  with  All  Cj  equal  to  zero.  With  all  *0,  L  in 
(2.1A)  is  zero  for  any  set  of  feasible  ;  hence,  L  is  at 
its  minimum  value  for  any  set  of  feasible  Xj .  Since  L  is 
at  its  minimum  value  for  any  feasible  set  of  x^ ,  any 
feasible  set  of  Xj  is,  by  the  above  definition,  a  solution 
to  the  linear  programming  problem. 

However,  we  must  not  only  find  a  feasible  solution  to 
the  linear  programming  problem,  we  must  also  find  a  positive 
feasible  solution  to  the  problem.  In  order  to  do  this,  we 
let 


Yi  y 


n+1 


J-l,2,...,n  (2.15) 


If  we  can  find  non-negative  values  of  yi>y2» • *  * ’^n+l 
which  satisfy 


n 

I  ^n+l^  "  ^i  . ^  (2.16) 

j-1 


then  Xj ,  as  defined  by  (2.15),  will  be  a  feasible  solution. 
If  we  can  somehow  assure  that  Is  positive,  then  all 

Xj  will  be  positive.  Rewriting  (2.16),  we  have 
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n 


I  ^j>'j 


i=l , 2 , . . . ,m  (2.17) 


If  we  now  specify  > ^2  ’  *  *  * ’^n+1  ’  have  a  linear  program¬ 
ming  problem  in  n+1  unknowns.  In  order  to  guarantee  that 
is  positive,  if  it  is  possible  for  it  to  be  positive, 
we  can  maximize  It  is  easy  to  see  that  we  can  maximize 

Yn+i  hy  setting 


I  -  -  (2.18) 

which  is  equivalent  to  setting  C, =C„=Ct=. . . =C  =0,  C  =  -1. 

If  the  solution  to  the  resulting  linear  programming  problem 
is  feasible  and  ^  0>  then  we  have,  by  (2.13),  a  positive 

feasible  solution  to  the  analogous  chemical  equilibrium 
problem  (1.1).  If  the  linear  programming  problem  is  feasible 
but  “  Oj  then  the  analogous  chemical  equilibrium  problem 

is  degenerate,  since  there  is  no  strictly  positive  solution 
to  the  problem.  However,  this  is  a  rather  trivial  kind  of 
degeneracy,  and  its  occurrence  usually  indicates  that  a 
mistake  was  made  in  setting  up  the  problem.  Hence,  this 
linear  programming  method  gives  us  a  way  of  finding  a  positive 
feasible  solution  to  the  chemical  equilibrium  problem  if 
the  chemical  equilibrium  problem  is  non-degenerate. 
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The  positive  feasible  solution  that  we  obtain  by  this 
method  will  generally  not  resemble  the  final  solution  of 
the  chemical  equilibrium  problem.  The  initial  positive 
feasible  solution  can  be  improved  by  the  following  tech¬ 
nique.  Define  b^^^  to  be  some  multiple,  between  zero  and 
one,  of  the  value  of  y  that  was  obtained  above.  Then, 
adjoin  to  the  linear  restraints  (2.17)  one  more  restraint 
of  the  form  *  ^m+l*  Ne.  t,  solve  the  linear  program¬ 

ming  problem  with  these  restraints  and  with 
...,  C  “c  ,  C  .  *0  (recall  that  the  lower-case  c's  here 
refer  to  the  c's  in  the  chemical  equilibrium  problem  (1.3)). 

The  solution  to  this  linear  programming  problem  will  give  a 
set  of  components  more  nearly  resembling  the  solution  to  the 
chemical  equilibrium  problem  than  did  the  components  calcu¬ 
lated  from  Eqs.  (2.17)  and  (2.18).  This  new  solution,  in 
turn,  may  be  improved  by  solving  another  linear  programming 
problem  (the  details  of  which  can  be  seen  in  SUBROUTINE  LP  in 
Appendix  A)  and  averaging  the  new  solution  with  the  old  solution 
In  order  to  solve  an  elaborate  chemical  equilibrium 
problem  it  is  not  sufficient  to  simply  use  a  method  which 
we  can  prove  converges  to  the  correct  solution.  Proofs 
of  convergence  generally  assume  infinite  computational 
accuracy,  but  since  we  are  usually  limited  in  practice  to 
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about  eight  significant  digits,  the  numerical  solution  will 
not  always  converge.  However,  it  has  been  observed  that 
the  closer  we  can  get  to  the  solution  by  the  initial  solu¬ 
tion  methods  described  above,  the  greater  will  be  the 
probability  that  the  numerical  procedure  will  converge. 
Furthermore,  not  only  will  the  probability  of  convergence 
be  greater,  but  the  number  of  iterations  to  get  to  the 
solution  will  be  fewer,  and  hence--when  an  improved  initial 
solution  is  used--the  computation  time  will  be  shorter. 
Unfortunately,  the  mathematical  methods  that  are  available 
for  analyzing  convergence  of  iterative  processes  do  not, 
in  the  case  of  the  chemical  equilibrium  problem,  enable  us 
to  prove  convergence  when  we  are  limited  to  finite  mathe¬ 
matical  accuracy.  Only  experience  with  a  particular  method 
will  tell  us  whether  it  is  a  useful  numerical  procedure 
to  use. 

In  the  next  section  we  consider  a  somewhat  more  general 
problem  than  the  chemical  equilibrium  problem.  This  prob¬ 
lem  is  considered  first  because  the  numerical  results  take 
on  an  especially  simple  form  when  the  additional  generality 


is  admitted. 
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3.  THE  LINEAR-LOGARITHMIC  PROGRAMMING  PROBLEM. 
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We  may  attempt  to  solve  this  problem  using  Lagrange 
* 

multipliers.  In  this  method  we  let 


and  then  set 


J 

for  j  =  1,2,3,...,N.  Performing  the  partial  differentia¬ 
tion,  we  get 


M 

c.  +  d.  log  X.  +  d.  -  y  Ti.a..  =  0  ,  (3.3) 

J  J  J  J  ^ 

i=l 

j=l,2,3, . . . ,N 


or,  when  rearranged. 


log  X.  =  d. 
J  J 


M 


y  n.a.  . 
L  I  ij 

i=l 


c  .  -  d  . 
J  J 


j=l,2,3,...,N 


(3.4) 


•k 

See  Kaplan, 
p.  140. 


Ref.-  3 ,  p. 


i28. 


or  Dantzif,  Rv^f.  2, 
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Exponentiating  both  sides  of  (3.4),  we. get 


X .  =  exp 
J 


M 


d.^  Vff.a  -  d.^c.  -  1 
3  L  1  i2  J  J 

j=l,2,3, . . . ,N 


(3.5) 


Note  that  for  (3.5)  to  be  a  solution  to  the  problem,  we 

must  have  all  x,  >  0.  We  assume,  in  the  remainder  of  this 

J 

section,  that  the  solution  does  have  all  x.  >  0.  Then, 

J 

the  problem  reduces  to  the  problem  of  determining  the  M  ir^ 
so  that  the  x^  from  (3.5)  satisfy  (3.2)  Equivalently, 
the  M  +  N  equations  (3.2)  and  (3.5)  must  be  satisfied  simul¬ 
taneously  by  the  proper  choice  of  the  M  +  N  unknowns, 

» ^2  ’  *  •  •  ’  consider  two  methods 

of  approximating  the  solution. 

In  the  first  method,  we  suppose  that  we  have  an  esti¬ 
mate  of  the  Xj  which  may  or  may  not  satisfy  (3.2).  We 
denote  this  estimate  by  y^ ,  and,  in  this  method,  solve 
Eqs.  (3.2)  and  (3.4)  simultaneously  by  making  a  linear 
approximation  to  log  x^ .  Since  we  have  the  estimate  that 
Xj  is  near  y^ ,  we  note  that  the  first-order  Taylor  ex¬ 
pansion  of  log  Xj  about  y^  is 
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x.-y . 

log  X.  =  log  y,  ^  ^  +  (higher-order  terms)  .  (3.6) 

J  J 


Dropping  the  higher-order  terms,  and  substituting  (3.6)  into 


(3.4)  and  solving  for  ,  we  have 


X.  =  y. 
J  J 


M 


y  TT.a.  .  -  d.^c.  -  log  y. 
J  A  1  i-J  J  J  ^  J 


-1 


i=l 


(3.7) 


j»l,2,3,. . . ,N 


Now,  if  we  substitute  these  x^  into  (3.2),  we  get 
M  /  N  \  N 

I  ( ■  "i  ^  "j  ^  • 

/-I  \j»l  ^  j-1 


i-l,2,3,...,M 


Denoting 


N 


il 


-  I 


,-l 

d .  a. .a  .y . 
J  i-J  12  3 


•(-■1 ,2,3,.  .  .  ,M 
i-1,2,3, . . . ,M 


(3.8) 


j-1 


and 


N 


®1  ■  ‘’l  yj  + 

j-1 

i-1,2,3, . . . ,M 


(3.9) 
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we  have 

M 

t-1,2,3,  .  .  .  ,M  (3.10) 

^“1 

Equation  (3.10)  is  a  set  of  simultaneous  equations  which 

can  be  solved  for  Tr-,T7  .  . 

i  z  n 

With  the  above  results,  we  can  now  define  the  iterative 
process  for  the  first  method.  At  each  iteration  we  have  a 
set  of  values  for  ,X2 ,  .  .  .  At  the  beginning  of  the 

iteration  these  values  are  called  , y2 » •  •  •  » Yjj*  the 

end  of  the  iteration  the  values  are  x^^  ,X2  , .  .  .  ,x^.  If 


is  small  for  each  j,  then  we  say  we  have  converged.  The 
magnitude  of  "small”  depends  on  the  nature  of  the  problem. 
If 


is  not  small  for  some  j,  then  we  have  not  converged  and 
the  iteration  must  be  repeated.  One  iteration  consists  of 
the  following  three  steps: 
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1)  Evaluate  terms  In  Eqs.  (3.8)  and  (3.9),  these 
terms  depending  on  , y2 »  •  •  •  » 

2)  Solve  Eq.  (3.10)  for  » •  •  •  » 

3)  Substitute  , .  .  .  into  (3.7)  to  get 

i  z  M 

For  this  problem,  in  this  generality,  we  can  say  noth¬ 
ing  about  whether  this  iterative  process  converges.  In 
the  next  section  we  will  show  that  the  chemical  equilibrium 
problem  is  a  special  case  of  this  problem,  and  one  for  which, 
with  appropriate  modification,  this  method  does  converge. 
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4.  THE  FIRST-ORDER  METHOD  FOR  SOLVING  THE 
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0 

if 

i 

m. 

j 

> 

n 

1 

if 

i 

> 

m. 

j 

n, 

and 

[j] 

■  i-m 

0 

if 

i 

> 

m. 

j 

n, 

and 

[j] 

^  i-m 

1 

if 

i 

> 

m. 

j 

> 

n, 

and 

i-m 

-  j-n 

0 

if 

i 

> 

m. 

j 

> 

n, 

and 

i-m 

^  j-n  . 

(4.4) 


For  all  j,  we  define 


d, 

J 


+1  If  j  ^  n 
-1  if  j  >  n  . 


(4.5) 


With  these  definitions,  it  has  been  shown  [4]  that  the  two 
problems  are  identical.  Next,  we  let 


X. 

J 


y .  +  0 . 
J  J 


(4.6) 


ffl  +  log  S.  +  1  . 
i  ®  i-m 


iim 


i>in 


Substituting  Eqs.  (4.1)  through  (4.6)  into  (3.7)  through 
(3.10)  and  simplifying,  we  have 
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I  “i/l  ■ 


’'j  ’fj]^ 


.  .  a  jf! 


(4.7) 


Za. .a  ,y, 

j-1 


ism 


I 


jc  <i-in) 


I  ^j>'j 

jc  <-C-m) 


i>in 


(4.8) 


■^>tn,  lim 


'>in,  i>in 


f  _ 


b.  +  y  a..y.(c.  +  leg  y.  -  1) 
1  L  ij^j'  J  ^ 

j-1 


^  yj(Cj  +  log  yy 
j  €  '' 


(4.9) 


y  ^  ^ ' 

L  I 


i=l,2,...,M  (4.10) 


The  directional  derivative  of  F  in  the  direction 
(0j^,02, .  .  .  ,0j^)  is  given  by  [1,  Theorem  8.11]  to  be 
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+  log^j)  . 

j=l 


N 


But,  if  we  compute  ^  ^  where  by  (3.7) 


j  =  l 


(4.11) 


^k+n  I  "m+k 


tr 

t 


-  log  S, 


-  1 


Su 

k  iTr+-k 


(4.12) 


k“l , 2 , . . . , p 


we  show,  in  Appendix  B,  that 


N  2 
OTd. 


n 


m 


n 


Z  ^  -  I  ^  “l/j  )  ■ 

j-1  J  j=l  l-l  \  J-1  / 


(4.13) 


Thus,  if  we  assume  that  (y^^ ,72 >  •  •  •  > Is  feasible,  we  get 

the  itijeresting  result  that  the  directional  derivative  of 

F  in  the  direction  (6, 9  )  is 

12  n 


j  =  l 


^  0 


j-1 


j 


(4.14) 


However,  it  is  also  shown  in  Api>  -)dix  B  that  the 
equality  on  the  right  side  of  (4.14)  holds  if  and  only  if 
the  values  for  y^  are  optimal.  We  further  note  that  if 
(y2^,y2.  •  •  •  »yn)  is  feasible,  then 
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n 

y  a.  .0 .  =  0 

j  =  l 

for  I  =  Hence,  if  ,  y2  >  •  •  •  >  is  feasible, 

then  (y-,+A.0.  ,  y-+A0_  ,  .  .  .  ,  y  +A0  )  will  be  feasible  for  any 
^11-^22  “^nn 

A  for  which  each  y.  +  A0.  is  positive. 

J  J 

We  now  state  the  first-order  chemical  equilibrium 
algorithm; 

1)  Calculate  (0^  ,©2 >  •  •  •  *6^^)  using  Eqs.  (4.7)  through 
(4.10). 

2)  Calculate  the  directional  derivative  of  F  in  the 

direction  (0^^  ,02  >  •  •  •  as  given  by  Eq.  (4.11); 

if  thir  quantity  is  not  negative,  we  are  done. 

3)  Calculate 

€  = 

f  is  a  number  that  represents  the  root-mean-square 
error  in  (y^^ , >2  >  •  •  •  > yj^)  •  ^  i-s  less  than  some 

given  numbv<jr  (say,  0.001),  we  are  done. 
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4)  Calculate  the  ratio  -y./0.  for  every  \  for  which 

J  J 

0.  <0.  Let  A,  be  the  minimuni  of  all  such  ratios 

J  1- 

and  let  A  =  min  (l,j3A^),  where  j3  is  a  number  less 

than  1  but  close  to  1  (say,  0.99).  We  now  per¬ 
form  the  following  steps  until  the  test  at  c)  be¬ 
low  is  satisfied: 

a)  Letz,=y.  -iA©.; 

J  J  J 

b)  Compute  the  directional  derivative  of  F  at 

z,  in  the  direction  (0-,0-,...,0  ):  f(A)  = 
j  1  2  n 

0  (c  +  log  z  ; 

J  J  J 

c)  If  f(A)  -s  0,  go  directly  to  step  5); 

d)  Replace  A  by  yX ,  where  0  y  <  1 ,  e.g.,  y  ~  ^ 

5)  Finally,  replace  by  y^  +  X0j  for  j  =  I,2,...,n, 
Steps  1-5  are  repeated  until  either  the  test  in  step  2  or 
the  test  in  step  3  is  satisfied. 

If  this  process  terminates,  the  solution  will  be 
optimal  within  the  specified  limits  of  accuracy.  It  may 
happen  that  the  process  does  not  terminate.  Since  the 

■k 

objective  function  F  is  convex  and  assuming  infinite 
computational  accuracy,  non- terminat ion  can  occur  only  be¬ 
cause  the  values  chosen  for  A  become  smaller  on  every 


Ref.  1,  Theorem  8.13;  Ref.  5. 
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iteration.  This  will  occur  only  if  some  is  approaching 
zero,  and  hence  (y^^ , y2  >  •  •  •  >  i-S  approaching  a  point  at 
which,  if  it  were  the  optimal  solution,  the  problem  would 
be  degenerate.  It  is  possible  for  this  to  happen  for  a 
non-degenerate  problem  for  which  the  initial  solution 
chosen  was  too  far  from  the  optimal  solution.  Cotivergence 
can  be  guaranteed  by  imposing  the  condition  that  the  value 
of  F  at  the  initial  solution  be  less  than  the  value  of  F 
at  any  feasible,  degenerate  point.  However,  it  is  not 
practical  to  impose  this  condition  on  the  initial  solution 
since  it  may  be  very  difficult  to  find  such  a  point.  In 
practice,  it  has  been  found  that  round-off  errors  cause 
more  difficulty  than  the  possible  selection  of  a  poor 


initial  solution. 
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5.  THE  LI NEIAR- LOGARITHMIC  PROGRAMMING  PROBLEJ^, 

SECOND-ORDER  METHOD 

In  the  first-order  method,  presented  in  Sec.  3,  the 

iterative  process  was  initiated  with  an  estimate  of  the 

value  of  ,  .  .  .  ,  In  the  second-order  method,  we 

assume  that  the  problem  is  as  defined  by  Eqs.  (3.1)  and 

(3.2),  but  that  we  have  initial  estimates  for  the  values 

of  ,  T7  .  .  ,  T7  .  Let’us  denote  these  estimates  by 
i  /  M 

A,,A  The  x.  can  then  be  evaluated  by  Eq .  (3.5), 

12  M  j 

substituting  A^  for  ir^.  These  x^ ,  however,  probably  will 
not  satisfy  Eq .  (3.2).  The  problem  of  the  second-order 
method  is  to  find  numbers  AA,  , , . , . , ,  such  that 

Tr.=A.-»-AA.  i  =  l,2,...,M  (5.1) 

111 


when  substiti ted  into  (3.5)  will  give  x^  that  satisfy  (3.2) 
In  order  to  accomplish  this,  we  first  use  the  x^ 
calculated  from  Eq .  (3.5)  to  get 


N 


=b.  -  ya..x. 

■i  1  ^  ij  J 


i  =  l,2, . .  .  ,M 


(5.2) 
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where  represents  the  amount  that  equation  i  is  in  error. 
Next,  we  evaluate 


bX 

I 


by 


dX 


r 


b. 

1 


N 

I 

j=l 


a . 


i-j 


exp 


dT^c. 

J  J 


N 


=  -  I 


,-l 

a. ,d .  x.a  . 

i-J  J  J  'O 


=  -  r 


li 


(5.3) 


j  =  l 


where  r^^  is  given  by  Eq.  (3.8) 
change,  dX^,  dX^,...,  in  Xj^,X2, 
is  given  by  dg^,dg2,...,  where 


^.=1 


If  we  make  a  very  small 
,  the  change  in  g^,g2, • 


i=l,2. 


-30- 


or 

M 

dgi  =  -  ^r.^dX^  .  i=l,2,...,M  (5.4) 

^=1 

We  would  want  dg^  to  be  equal  to  -g^  as  computed  by 
Eq .  (5.2).  If  we  make  the  approximation  that 


ax . 

j 


is  constant  over  the  domain  considered,  we  can  set 
~  ^ I '>  write 

M 

gi  =  •  i=l,2,...,M  (5.5) 

l=l 

Equation  (5.5)  consists  of  M  equations  in  the  M  unknowns 

AX^  ,  ^2  ,  .  .  .  ,  We  may  thus  solve  Eq  .  (5.5)  for 

AX,  ,  AX„  ,  .  .  .  ,  AX„  and  compute  n.  ,rr.  ,  ,  .  .  ,  n..  from  (5.1).  If 
1.  z  M  1  Z  M 

the  assumption  about 


being  constant  over  the  domain  considered  was  correct,  then 
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the  computed  from  (3.5)  with  these  values  for  ir^  will 
satisfy  (3.2).  However,  in  general,  they  will  not  satisfy 
(3.2),  but,  if  we  were  close  enough  to  the  solution  so 
that  the 

!!i 

dX 

did  not  vary  greatly  in  the  domain  considered,  then  the  new 

values  for  should  come  closer  to  satisfying  (3.2)  than 

did  the  first  set  of  x.. 

J 

With  this  assumption,  we  may  now  state  the  iterative 
process : 

a)  Using  the  values  at  hand  for  , .  .  .  ,  Ww, 

I  Z  n 

evaluate  (3.5). 

b)  Using  the  values  for  x^  obtained  in  step  a^ 
evaluate  (5.2).  If  the  jg^|  are  sufficiently 
small,  we  are  done. 

c)  Compute  r.  using  (3.8)  and  solve  (5.5)  for  AX.. 

d)  Denoting  the  it ^  in  step  a  by  A^,  we  get  new 
by  (5.1). 

Steps  a-d  are  repeated  until  the  |g^l,  computed  in  step 
b,  are  sufficiently  small,  or  until  they  show  no  more 
improvement . 
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There  is  no  proof  of  convergence  for  this  method. 

In  fact,  the  method  presented  here  is  unlikely  to  converge 
unless  the  starting  values  of  »  •  •  •  » very  good, 

and  even  then  there  may  be  no  convergence.  This  method  may 
be  used  on  the  chemical  equilibrium  problem  after  the  first- 
order  method  has  resulted  in  a  reasonably  good  solution. 

If  the  obtained  from  (3.10)  in  the  final  iteration  of 
the  first-order  method  are  used  to  initiate  the  second-order 
method,  the  av.curacy  produced  by  the  second-order  method 
will  generally  be  better  than  that  which  could  be  achieved 
by  use  of  the  first-order  method  only. 
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6.  THE  SECOND-ORDER  CHE>1ICAL  EQUILIBRIUM  ALGORITHM 

In  ord^r  that  the  second-order  linear-logarithmic 

method  be  set  in  the  form  of  a  chemical  equilibrium  problem, 

the  same  definitions  as  given  in  Sec.  4--i.e.,  Eqs.  (4.1) 

through  (4.5)--are  used  here.  Since  the  second-order  method 

is  best  used  after  the  first-order  method  has  been  applied, 

the  initial  values  of  ir .  for  the  second-order  method  must 

1 

be  specified.  The  first-order  method  gives  a  set  of 
which  are  related  to  rr^  by  Eq.  (4.6).  The  computed  by 
means  of  (4.6)  are  appropriate  initial  values  for  the  second 

order  method.  Using  these  initial  values  for  n^,  the  second 

order  chemical  equilibrium  algorithm  is  an  iterative  process 
for  which  each  iteration  consists  of  the  following  steps: 

1)  Using  the  current  values  for  (tt,  , ,  .  .  .  , ff„)  , 

1  z  M 

evaluate  ,  .  .  .  ,  by  means  of  (3.5). 

2)  Calculate  ,g2  >  •  •  •  >  t>y  means  of  (5.2)  and  set 

3)  Compute  r^^  from  (4.8)  and  solve  (5.5)  for 

AX  )  AX  2  ,  .  .  .  ,  ■ 

4)  Let 

M 

P  =  max  '  AX  .  '  . 

.  ,  I 
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If  P  <  6,  where  6  is  a  small  positive  number  such 
as  10  we  are  done;  otherwise,  let  Q  =  min 
5)  Replace  by  -f  Q  for  i  =  1,2,...,M. 

Steps  1-5  are  repeated  until  the  test  at  4)  is  satisfied. 

P  should  decrease  at  every  iteration;  however,  when  the 
values  for  get  close  to  their  optimal  values,  P  may 
not  become  zero  due  to  round-off  error.  In  order  to  prevent 
an  endless  repetition  .of  steps  1-5  due  to  the  selection  of 
too  small  a  6 ,  we  can  test  P  against  the  value  of  P  at  the 
previous  iteration.  If  this  value  has  increased  over  the 
previous  iteration,  it  can  be  assumed  that  this  method  has 
obtained  as  accurate  a  solution  as  possible,  and  we  can 
terminate  the  iteration  process.  The  reason  for  inserting 
the  factor  Q  above  is  to  prevent  the  from  varying  too 
much  on  one  iteration. 
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7.  SUMMARY  OF  THE  COMPUTATION  PROCEDURE 

The  best  method  for  starting  the  solution  of  the 
chemical  equilibrium  problem  depends  on  whether  an  estimate 
for  the  solution  vector  is  available.  The  projection  method 
should  be  used  when  the  problem  being  solved  is  a  slight 
variation  from  a  problem  previously  solved,  and  in  this 
case,  the  values  used  for  in  (2.9  -  2.12)  should  be  the 
solution  vector  to  the  previous  problem.  Even  when  the 
estimate  is  no  better  than  an  intuitive  guess,  the  pro¬ 
jection  method  may  still  be  used.  The  linear  programming 
method,  then,  may  be  used  as  a  back-up  if  the  projection 
method  produces  a  non-positive  component.  Of  course,  if 
no  estimate  is  available,  the  linear  programming  method 
would  be  used  immediately  to  provide  an  estimate. 

The  recommended  procedure  is,  then,  to  use  the  first- 
order  method  until  either  no  further  progress  can  be  made 
with  this  method  or  until  the  amount  of  change  becomes 
small  from  iteration  to  iteration,  and  then  to  use  the 
second-order  method.  It  has  been  found  that,  for  reason¬ 
ably  large  problems  (say  m  =  30,  n  =  100),  the  point  at 
which  progress  ceases  in  the  first-order  method  usually 
occurs  when  the  indicated  corrections  to  the  components 
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of  the  solution  vector  average  about  one  per  cent  of  the 
components;  that  is,  when  (3.5)  is  accurate  to  about  two 
significant  digits.  A  switch  to  the  second-order  method 
at  this  point  usually  yields  quite  accurate  results  in  two 
iterations  of  the  second-order  method.  The  second-order 
method  usually  satisfies  (1.1)  to  an  accuracy  of  about 
five  significant  digits  on  a  machine  that  carries  eight 
significant  digits.  This  accuracy  is  typically  about  three 
orders  of  magnitude  above  what  is  usually  obtained  in 
experimental  data. 

To  summarize,  the  typical  procedure  for  solving  a 
chemical  equilibrium  problem  is  the  following: 

1)  If  an  estimate  is  available,  use  the  projection 
method  to  obtain  a  feasible  estimate. 

2)  If  step  1  yields  a  strictly  positive  estimate,  go 
to  step  3,  but  if  tt.e  projection  method  yields  non-positive 
components,  or  if  there  was  no  initial  estimate,  then  use 
the  linear  programming  method  to  get  an  estimate. 

3)  Use  the  first-order  method  until  one  of  the  tests 
described  in  Section  4  is  satisfied. 

4)  Use  the  second-order  method  as  described  in  Section 


6. 
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Appendix  A 

A  FORTRAN- IV  PROGRAM  FOR  SOLVING  THE 
CHEMICAL  EQUILIBRIUM  PROBLEM 

GENERAL  DESCRIPTION 

The  program  described  here  is  a  set  of  FORTRAN-IV 
subroutines  for  solving  chemical  equilibrium  problems. 
The  calling  sequence  used  is  merely  the  statement: 

CALL  SOLVE 


Communication  of  data  into  and  out  of  the  subroutines 
is  accomplished  by  a  block  common  statement: 

COMMON/ SLVE/ IV (30) ,T0L(20) , NR(55 , 2) , B(55) ,KN(120) ,X(121) ,C(121) , 

1  KL(26),NAM(25,2),A(55,121),PIE(65),V1(65),V2(65),V3(65), 

2  V4(65) ,XMF(120) ,X1(121) ,X2(121) ,X3(121) ,XBAR(25) , R(65 , 65) 

The  data  that  must  be  input  before  CALL  SOLVE  is 
executed  consist  of  the  following: 


COMMON  Location 
IV(1) 

IV(2) 

IV(3) 

IV(4) 

IV(6) 


Quantity 

m 

M  (  =  m+p) 

P 

n 


Number  of  the  output  unit. 
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COMMON  Location 
IV(7) 


IV(9) 


Quantity 

Print  flag:  -1  =  minimal  amount  of 
messages;  0  =  one  message  per  itera¬ 
tion  step;  +1  =  all  messages. 

Maximum  number  of  iterations  to  be 
allowed . 


B(i) 

X(j) 


b^j  1.  “  m • 

y.,  j=l,2,...,m,  where  y,  is  the 
Initial  estimate  of  the  solution. 
If  no  estimate  is  available,  set 
X(J)  =  0. 


C(j) 

A(i,  j) 

In  addition, 


c.,  j— ‘l,2,«*.,n. 

j 

a^j,  L^l,2,...,m;  j~l,2,«.«,n. 
all  components  in  one  compartment  must 


have  consecutive  subscripts.  That  is,  components  1,2,3,.. 
t  be  in  compartment  1;  components  k^+1 ,  kj^+2,  ...,  k2 


mus 


must  be  in  compartment  2;  ...;  and  components  k^  » 
kp  f  •••>  ni'JSt  be  in  compartment  p.  These  k's  are 
communicated  to  the  subroutines  by  setting 


KL(1)  =  1 
KL(2)  =  k^+l 
KL(3)  =  k^+l 


KL.(p)  =  kp_^-Pl 
KL(p+l)  =  kp+1 
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In  other  words,  KL(k)  is  the  number  of  the  first  component 
in  compartment  k,  and  KL(p+l)  is  equal  to  n+1. 

The  above  are  the  only  numbers  that  need  be  set  in 
order  that  CALL  SOLVE  will  solve  the  chemical  equilibrium 
prv^bLm.  However,  in  order  that  the  program  can  write 
messages,  in  cases  of  infeasibility,  etc.,  names  for  the 
rows,  components,  and  compartments  may  be  input: 


COMMON  Location 
NR(I,1),  NR(I,2) 
KN(J) 


NAM  (K,l),  NAM(K,2) 


Quantity 

Two-word  row  name  for  row  I, 

One-word  component  name  for 
component  J. 

Two-word  compartment  name  for 
compartment  K. 


In  addition,  TOL(l)  through  T0L(5)  are  tolerances  used 
by  the  program.  If  they  are  zero  when  the  program  is 
entered,  they  are  set  by  the  subroutines  to  nominal  values. 
These  values  may  also  be  set  by  the  user  of  the  subroutines 
in  which  case  the  nominal  values  will  not  be  set  in  the  sub 
routines.  These  tolerances  are  the  following: 

Nominal 

Tolerance  Value  Meaning 

c  in  step  3  of  the  first- 
order  method  (see  Sec.  4). 


TOL(l) 


0.01 
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Tolerance 

Nominal 

Value 

Meaning 

T0L(2) 

10“^ 

6  in  step  4  of  the  second- 

order  method  (see  Sec.  6) 

T0L(3) 

Minimum  value  any  x.  is 

J 

allowed  to  have. 

T0L(4) 

10-^ 

Minimum  starting  value  that 

any  component  will  have  is 
the  lesser  of  T0L(4)  and 


^^n+1  ■ 

-  8 

T0L(5)  10  Problem  is  assumed  to  be 

degenerate  if  any 
becomes  less  than  T0L(5). 

With  the  above  as  input,  the  statement  CALL  SOLVE  will 
cause  an  attempt  to  solve  the  chemical  equilibrium  problem. 
If,  upon  completion  of  this  attempt,  a  solution  is  obtained, 
the  cell 

IV(IO) 

will  contain  a  1  and  the  following  data  will  be  in  storage: 

Data 

x^,  i=l,2,.,.,n  (the  solution). 
k=l,2, . . . ,p. 

,  i  =  l ,  2  ,  .  .  .  ,m. 

X.,  i^l,2,...,n. 


COMMON  Location 
X(i) 

XBAR(k) 

PIE(i) 

XMF(i) 
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If  IV(IO)  is  not  1,  the  subroutines  have  failed  to  solve 
the  chemical  equilibrium  problem.  The  reason  for  this 
failure  is  written  on  output  unit  IV(6) .  In  such  a  case, 

X(i)  will  contain  the  latest  value  of  these  quantities. 

SUBROUTINES 

There  are  nine  subroutines  in  the  set  used  for  the 
solution  of  the  chemical  equilibrium  problem.  A  brief 
description  of  these  subroutines  follows. 

1.  Subroutine  SOLVE 

SOLVE  is  the  master  subroutine,  and  is  divided  into 
four  functional  segments.  Each  segment  calls  other  sub¬ 
routines  which  do  specific  tasks.  The  four  segments 
are: 

a)  The  projection  and  linear  programming  routines 
for  obtaining  the  initial  solution  (lines  18-42) . 

b)  The  first-order  method  (lines  43-122). 

c)  The  second-order  method  (lines  123-163). 

d)  Output  messages  (lines  164-203). 

2.  Subroutine  BAR 

BAR  calculates  the  S,  . 

k 
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3.  Subroutine  BERRQR 
BERROR  calculates 

N 

g.  =b.  -  ya..x.  .  i=l , 2 , . . . ,M 
1  L  Z.  ij  J 

j  =  l 

4.  Subroutine  DEL 
DEI  sets 

m 

w.=  )a.,q,  .  1=1,2,., .,n 

j  L  ij^i 

i  =  l 

5.  Subroutine  RCALC 

RCALC  calculates  the  array  (4.8). 

6.  Subroutine  CLOG 
CLOG  computes 

A 

a.  =  c.  +  log  X.  .  1  =  1,2,. ..,n 

J  J  J  j  >  >  > 

7.  Subroutine  LP 

LP  sets  up  the  linear  programming  problems. 

8.  Subroutine  SIMPLE 

SIMPLE  solves  the  linear  programming  problems. 


Information  is  communicated  to  this  routine  via  a 
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calling  sequence  rather  than  by  COMMON  as  in  sub¬ 
routines  1-7.  The  dimension  of  A  in  SIMPLE  should 
agree  with  the  dimension  of  A  in  the  first  seven 
subroutines,  but  all  other  dimensions  are  dummy 
statements . 

9.  Subroutine  MATINV 

MATINV  solves  simultaneous  equations.  As  in 
SIMPLE,  no  COMMON  is  used.  The  dimension  of  A  in 
MATINV  should  agree  with  that  of  R  (not  A)  in  SOLVE. 
All  other  dimensions  are  singly  subscripted  and  are 
irrelevant  as  to  magnitude. 

*  *  * 

Each  of  the  first  seven  subroutines  has  a  COMMON 
statement  which  should  be  the  same  in  all  seven.  The 
dimensions  of  the  variables  in  this  COMMON  statement  may 
be  set  to  the  values  for  the  largest  problem  to  be  solved. 
With  m,  M,  p,  and  n  as  previously  defined,  these  dimen¬ 


sions  must  be  at  least: 
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Symbol 

Minimum  Dimens 

IV 

30 

TOL 

20 

NR 

(m,  2) 

B 

m 

KN 

n 

X 

n+1 

C 

n+1 

KL 

p+1 

NAM 

(P,2) 

A 

(m, n+1) 

PIE 

M 

V1,V2,V3,V4 

M 

XMF 

n 

X1,X2,X3 

n+1 

XBAR 

P 

R 

(M,M)  . 

A  listing  of  these  subroutines  follows.  This  listing  does 
not  necessarily  represent  an  actual  program.  The  language 
used  was  that  version  of  FORTRAN  described  in  [6'.  The 
machine  used  for  the  solution  of  chemical  equilibrium 
problems  was  the  IBM-7044,  which  uses  a  floating-point 
number  with  eight  bits  for  the  exponent  and  28  bits  for 
the  sign  and  mantissa. 
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LISTING 


subroutine  solve  souoi 

COMMON/SLVE  /  IV(30)tTOL(2u)tNR(S5.2)#B(S5)t<.N(12wi).X(12nfC(121lt  SCC02 

1  KH  26) #NAM( 25  t2 )tA(55»121)»PIt<65).Vl(65)tV2(65)tV3t65)t  S000  3 

2  V4( 65) ,XMF I120)tXl(121)tX2(121).x3(121) tXbAR (25),Rl65t65)  S0004 

INTEGER  PF  SC005 

equivalence  (TOL! 3).XMlN)t(TOLlA).XSTART).(TCL(5)  ,bARf/IN )  S0006 

equivalence  (  I V( 1  ) tM)  , (  I  V ( 2 ) .menu ) t ( I V( 3 )  .NCQVP )t(IV(4).NtNT0T),  oJCOT 

1  (IV(5).NlT)*(IV(6).NOT)t<IV(7),PF),(lv(8).ITER).(IV(V).lTi/AX).  S0C06 

2  ( I  V  I  )  t lERROR ) . (  IV(  1 1 )  .LASTCP )  .  (  I  VI  12 )  .XE )  S0039 

DIMENSION  DX  (  1  )  .ALPHA  (  1  ),  Trt  (  1  )  ,(j(  1  )  SOOIC 

equivalence  ( G. VI )  . < DX ,X 1 )  .  ( AlPhA  .X2  )  .  ( TH.X3 )  SOOll 

IF  ( TOL( 1 ) .LE.0.0 )  TOL(l)  =  C.Cl  S0C12 

IF  ( TQL ( 2 ) .LE  .C.O )  TOL(2)  =  l.E-5  S0013 

IF  (XMIN.LE.-.O  )  xy.lN  =  l.t.-12  SC014 

IF  (bARMIN.LE.C.O)  cARKIN  *  l.E-0  ouG15 

IF  (  I  TMAX  .  Lt. .  -  )  ITFAX  *  4u  S0016 

DO  152  J  =  1.  NTuT  S0C17 

IF(X(J).l£.>..)GQTC5  sc  018 

152  CONTINUE  oS019 

C  IF  X  IS  STRICTLY  POSITIVE.  oEGlN  PROJECTION  SCC20 

CALL  bAR(  X.XbAR  )  S0C21 

2  CALL  dERRCR(ERR)  SDC22 

call  RCALC  S0023 

call  mat  I  NV  (  r' .VENu.G.-l  .  V2  t  V3  •  V4  .KE  )  S0024 

IF  (XE.NE.v)  GO  TO  5  SCC25 

CALL  del  (DX.O)  SGC26 

DO  3  X  =  1  .i.CuMP  SG027 

XTA  =  xl(X)  S032a 

XTB  =  XL ( <♦!  )-l  oG029 

MX  =  M  X  S0030 

Do  4  J  =  XTA.XTb  SCC31 

X(J)  »  X(J1  *  (  1.  ♦  OXIJ)  ♦  G(VX)  )  SGC32 

IF  (X(J).lE.u.)  go  to  5  SGC33 

4  continue  S0034 

3  CONTINS-  3DC35 

GO  TO  7  S0036 

C  LINEAR  PROGRAMMING  RCuTInE  SCC37 

5  call  LP(XF,  S0C38 

IF  (XF.NE.D)  GO  TO  10oC6  S0039 

7  call  OARIX.XoAR)  S004C 

CALL  CLOG(X.xuAR)  jCCAI 

FE2  »  l.E>2u  S0C42 

C  first  order  METhQO  loop  S0043 

DO  899  ITER=1.ITMAx  S0CA4 

call  BERROR(ERR)  SCC45 

DC  711v  1*1. mend  Su046 

PIE(I)  *  So047 

7110  continue  S004fa 

DO  7111  X  *  1.  NCO.MP  jvc49 

XTA  *  XL(X)  SDC5C 

XTa  =  XL(X>1 )  -  1  ouC51 

MX  =  M  ♦  X  S0^52 

Du  7112  J  *  xTA.  XT3  o  53 

AX  «  AlPhA(J)  *  X(J)  SDC54 

PIE  (MX)  *  PIE(N.X)  ♦  AX  SOD  5  5 

DO  7113  I  *  1 .M  S0056 

PIE(I)  *  PIL(I)  ♦AX  *  A(i.J)  S0057 

7113  CONTINUE  30058 

7112  CONTINUE  SCC59 

7111  CONTINUE  Sj060 
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DO  7114  I  «  l.,1 

S0C61 

PIE( I  1  »  G( n  ♦  PIE( I  1  • 

SOC  62 

7114 

continue 

o0063 

call  rcalc 

S0C64 

CALL  MAT  INV  (S  tVESi^.PI  E  ,-l  ,V2  •  V3  tV4,<E  1 

SCC65 

IE  ( ICE.NE  »  GO  TO  K003 

SC066 

DMAX  «  l.E*2C 

O-067 

call  0£L(Tri,Pl£i 

SCC68 

71^5 

GNORM*w* 

iOC69 

TOA  «  0. 

S007C 

EE  =  0. 

S0C71 

DO  71^4  iC«ltNCOMP 

3CC72 

MA  »  M  ♦  < 

SCC73 

KTA  s  KL 1 <) 

S0C74 

kTo  *  <L(A»1)  -1 

S-075 

DO  7103  J  «  KTA,  k,T3 

ju076 

Th(J)  *  Trt(J)  ♦PI£(M<|-  ALPhA(J) 

S0077 

GNORM  *  GNOrM  ♦  Tm(J)  **2 

S0078 

TH( J)  a  TH( J)  *  X( J) 

00079 

TOA  a  TOA  ♦  Tm(ji  *  ALPriAIJI 

-CC  80 

IE  (  X  (  J  I  ,LT  .-OMAX*TH(' J  n  UMAX  =  -X(J)/IHIJ) 

SOCBl 

EE  a  F£  ♦  X(J)  *  ALPrtA I J ) 

oOC82 

7103 

CONTINUE 

S0C83 

7104 

continue 

SG0d4 

EPS*  SORT  (  GTOR^/ElOAT  (NT^T)  ) 

;>008  5 

DEE  a  FE  -  FE2 

-  jC  66 

EE2  *  EE 

8- JC  7 

IF  (I  TER.t-u.l  )  GO  TO  712- 

icoaa 

ITR  a  iter  -  1 

8CC89 

IFIPF.GE.-)  wRI TE (NOT.799  )  ITR,  0E£,CPTL.£PS 

SCG90 

712v 

CPTL  aAMlf.l  (  1.,  ,99*-<MAx  ) 

SCU91 

IE(PF.GT.j)«RITr (NCT,d24l)  D^A a ,0PTL , T JA ,ERK 

8jC92 

IF  ( EPS.LE.TOL ( 1 ) )  GO  TO  0^69 

Su093 

-li-i 

IE  (TOA,o£,m,)  go  to  6a.« 

o0094 

i\ »  1 

DO  8265  II  *1,54 

S0C95 

DO  63-1  J  *  1,N 

8-096 

OXtJ)  *  AMAX1(X(J)  ♦  LPTL»rn(J)  ,XVIM 

8-097 

8301 

continue 

30098 

call  L>AR(t/X,X  — AK) 

30099 

call  Cl-u ( ox  ,XuAR ) 

30100 

TOA  a  ^ , 

30101 

DO  8266  J  *  1,NT0T 

30102 

TOA  a  TOA  ♦  Tn( J I *ALPnA ( J 1 

30103 

8266 

CoNTINUt 

S0104 

IE  (PE  .GT.  -  1 /.R  1  T_  (  NOT  ,  8  262  )  I  I  ,C?T  L  ,  TOA 

30105 

IE  (  TDA,LT,0,)  go  to  828 

301  06 

8264 

OPTL  *  -PTL  /I. 4142 

80107 

8266 

continue 

30108 

call  8AR(X,XbAR) 

30109 

GO  TO  8271 

30110 

828 

DO  8281  J  *1,NT-T 

30111 

X( J)  a  0X( J) 

30112 

8281 

continue 

30113 

EE  =  0. 

30114 

00  8231  J*1,N 

30115 

EE  a  FE  ♦  alpha ( J ) ax ( J ) 

30116 

8231 

continue 

30117 

8288 

call  SSAiTCH(5,LAi£L) 

30118 

IE  (LAdEL,NC.2)  GO  TO  100-4 

30119 

bV9 

CONTINUE 

30120 

47 


C  EiNO  OF  FIKST  Gi<DER  METHOO  LOOP  S0121 

GO  TO  00122 

600O  ITERl  »  ITER  ♦  1  a0123 

PMAX  *  l.E+2u  o0124 

PMAXl  =  l.E-*-21  30125 

C  SECOND  ORDER  y.ETHOD  LOOP  S0126 

DO  60^2  ITER  *  ITERltITv-AX  S0127 

CALL  DELlDXtPIE)  S012a 

DO  6u.3  x  =1.NC3"P  S0129 

MTA  =  XL(K)  S0130 

PTB  X  <L(<^1)  -  1  S0131 

DO  6wl-  J  X  PTAt'^Te  S0132 

XMF(J)  xExP(DX(J)-CIJl)  SO  133 

X(J)  X  xf'F  (  J  )  *XaA.R  (  <1  SC134 

601vj  CONTINUl  S0135 

IF  ( XtJAR  (  X )  ,LE.  jAR-1  IN  I  00  TO  l-oC5  S0136 

Sv-v/S  continue  S0137 

IF  (  P-1AX  .Lc.  .TOL  I  2  )  .OR.  1  p  lAX.Gt  .P*''AX1  .ANu.P:'AX.GE.PKAx2  )  )  S013U 

1  GO  TO  1>^^--1  o0139 

call  SERROk(ERR)  SC140 

60v.6  CALL  RCAlC  S0141 

call  XAT INVIR.vENOtG.  - 1 . V2 t V3 . V4 ,K £ )  S0142 

IFIXE.NE.^)  GO  TO  1-003  SG143 

PMAX2  =  P1AX1  S0144 

PP.AXl  X  PMmX  i0145 

PPAX  X  SOI  46 

DO  6^04  I  =  IfPENJ  S0147 

P.VAX  xAMAXl  (  P-AX,  Aio  (G(I))  )  30146 

6Jv4  continue  S0149 

IF  (  PMAX  .  E  0 .0  .  -  )  GO  TO  1><-01  30150 

ZM  xAMINl  (  l./PMAX.l.)  S0151 

DO  6-i05  I  =ltM  S0152 

PIE(I)  »  PIE(I)  ♦  ZM*  Oil)  30153 

60v-5  continue  30154 

DC  6-1 1  K  x  1  ,f4C0PP  30155 

MX  X  K^x  30156 

XOAR(X)  *  XbARlX)*  EXP  (  2 X  *  GlMX)  )  30157 

6-11  continue  30156 

IF  (PF.GE.-)  wRI TLl NOT .6-99 )  I T £R .PMAX .ERR  30159 

call  33wTCm( 5 .lAuEL )  SO160 

IF  lLAdEL.NE.2)  GO  TO  1-004  S0161 

60-2  CONTINUE  30162 

C  END  OF  3t.C0ND  ORDER  VETmOD  LOOP  30163 

1- --2  lERROR  X  2  30164 

WR I TE ( NOT .2-D-2 )  S0165 

20--2  FORMAT(27rt  ITERAION  LIMIT  EXCEEDED  )  30166 

ITER  *  ITMAX  30167 

GO  TO  1-000  S0168 

10003  lERROR  x  3  -0169 

•RITE(NOT.2---3)  X£  S017C 

2- 0-3  FORMAT(21M  R  MATRIX  HAS  NULLlTy.I3)  o0171 

GO  TO  1-Oj-  30172 

1- --4  IERROR  *  4  5-173 

V.RITEIN0T.2---4)  S0174 

2- --4  F0KMAT(56M  SOLVE  ROUTINE  TER'  INATE-  DECAu'SL  SENSE  SSlICrl  5  IS  DO/.N  30175 

1)  30176 

GO  TO  100-0  3C177 

10005  IERROR  x  5  30176 

WRITEINOT. 2-0-5)  NAM ( X . 1 ) . NAM ( X  .  2  )  30179 

2-0-5  F0RMAT(13M  compartment  .2A6.10M  TOO  SMALL  )  30180 
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LASTCP  *  < 

GO  TO  KOJO 

lERKOi^  =  6 
GO  TO  1  J  V  w 
1  Ouw  1  1  EWSOi'  =  1 

l^Owi-  RETURN 

8i;4l  F0RMAT(13M  vax=  IPEI  2  *4  t  1  3M,  opt  LA’^ODAx  i,  1C.  3  *6^  ♦  TDAxt.12 

1.5.16^.  MAX  Rwa  EARwRxE12.b) 

8267  IF  '3F.oE.-l  aRUE  (NOT, 8268)  ITLR 

8260  CR.VAKKH  I  T  E  R  AT  I  ON  .  I  4  .  i  POSIIIvE  ICA,  GO  lO  veimOD  2  ) 

GO  TO  6C 

6cb't  IF  (PF.GE.-)  aRITi.  (NOT. 827  )  ITER 

b27v-  FoRF.AT  (  1-H  1  TlRAT  I  ON.  I4.42n  AV  TmlTA  LE68  I  nAN  I  JL  1  1  )  •  Gu  lu  KEThc 
10  2) 

GO  TO  bC:-, 

8271  IF  (PF.GL.^)  aRITl  (NOT, 6272)  ITER 

8272  FuRF'AT(1Cm  I  TERAT  1  on,  1  ,  36n  jTEP  OlZE  TOO  GMALL  »  GO  Tij  VETHOD  2) 

GO  TO  6000 

8262  Fv.RmaT(1-X,  4mSTEP,I2,  9h  L  A  ■•'uO  A  x  i  PE  1  C  .  3  ,6h  ,  TDAxtlG.b) 

797  Fu«MAT(lv«  ITERATION. I4,24M  Cm«NoE  IN  FRL-  ENE RG Y x i PE  1 9 , b • 1 2H 
ISTEP  GI2E-E16.8.1  H  AV  ThE  T  Ax  E  1  2 . 5>  ) 

6J99  FCRmaT(Uh  I  TERAT  I  on.  1  4. 1  9n  ..'AX  CHANGE  IN  P  I  E  =  1  PE  1  G  .  8 . 1  5h  vax  KOa 
1  ERROR  =  E 1 3 . 8 
END 


60183 

30184 

60183 

60186 

30187 

30188 
60189 
50190 

30191 

30192 
6ul93 
30194 

30193 

30196 

30197 

30198 

30199 
602  JO 

30201 

30202 

30203 

30204 
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subrout  I  \E  bAR  ( ,  v.3AR  )  'rtOCOl 

COMMON/SLVE/ I V( 33)»TOL(2-)*NR<55t2)»B(S5J.KNI12U),X(121).C(121)»  *U002 

1  KL(26)»NAM(25»2)*A(55»12n»PIE(65)»Vl(65)»v2i65)iv3i65)*  V*UUU3 

2  V4(63)tXMF(120)»Xl(121).X2n21).X3(121)  »XBAR(  25  )  .R(  65.65)  ^ftCOOA 

equivalence  (IV(1).M).(IV(2)  .MEND  )  .  I  1  V  (  3  )  .NCOMP  ).(1V(4).N.M0I).  ^/CC05 

1  (lV(5).NIT).(IV(6).NOT).(IV(7).PF).(lV(8).IIcR).(IV)9).ITMAX).  .VC  006 

2  ( I V( 10) . lERROR  )  .  (  1V( 11 ) .LASTCP ) . ( 1V( 12) .K£)  WU007 

DIMENSION  W(1).WBAR(1)  W0006 

7  DO  701  K  «  l.NCOMP  *0009 

<TA  »  KLIX)  v.OClO 

KTB  «  <L(K+1)  -  1  UOOll 

*'BAR(K)»0»  *0012 

DO  702  J  =  XTA.KTB  *0013 

'*bAR(X)  =  WBAR(<)  ♦W(J)  *0014 

7j2  continue  WQ015 

701  CONTINUE  *CC16 

END  .VOOIT 
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SUoKO^'TI'JE  (  -  'AX  )  bOCOl 

CON'-;v^.\/,.LVE  /  r.M  3^  I  .  TOL  (  2  J  htN'^  (  55 . 2  )  .  M  ';'■>)  .fCiMl  1  i;u  )  .X  I  1  2  1  )  *C  I  12  1  )  t  bU'J'J2 

1  nL  (  26  )  ,,\A,v,  (  2  5.2).A(55.121).Plt(65).yl(65).V2(o5)tV3l65).  bCC03 

2  V4(  65  I  ,X---F  (  12  .  )  *  XI  (  121  )  .X2  (  121  )  .X^  (1  2  1  )  .x‘,  (  25  )  ,«  (  65 ,65  )  bOnOA 

Ei^UlVALCMC:  (IV(n»''l,(lV(21,  -lEN  j  )  ,  (  iv  (  3  )  )  ,  (  I  v(  4  )  ,N  ,MOT  )  ,  BC^CCS 

1  (IV(5),MIT),(IV(6)»iSwT|,(IV(7).PF),ilviei*liiS),iIvi9),Ii  i-.Aa  )  ,  bUv06 

2  (  I  V  (  12  )  ,  I  Liiiric.9  )  ,  (  U  t  1 1  )  ,L^  J  t  CP  )  ,  l  I  V(  12  )  »'■;:  )  dUC3  7 

DIi'iLNSlO.'J  G(  1  )  BwCOa 

EuUIVAlL'JCc  (o,V1)  d0C09 

D01^1I  =  l,''l  oUOlO 

jGCll 

Uulv^2J=l,i'l  duC12 

IF  (A(  I  »J)  )  2T  =  2T  -  X(J)  *  A(1»J)  bC013 

1'>2CCNTIiN1JE  bOClA 

G(  I  )  =  ZT  b(  1  )  bJ015 

U1  COMIiMUE  bCC16 

GO  11-.  K  =  l.NCOV,?  bo017 

2T  =  V,.  uGCia 

MTA  =  <L(x.)  b3C19 

MTB  =  <L(<>1)  -  1  dUCZC 

DO  111  j  =  ;-;ta,:-,T6  60021 

ZT  =  ZT  +  X(J|  uOCZ2 

111  continue  BU023 

Mk  =  ,V|  s  b0024 

G(nN)  »  XuAR(i;)  -  ZT  BOO 2 5 

llw  CONTINUE  60026 

BMAX  =  60027 

DO  120  I  .  l.MEND  60026 

IF  (  ABS  (  G  <1  )  )  .GT  ,  AoSIB'IAXj  )  2'-AX  r  G(I)  B0029 

12'^  CONTINUE  6o030 

RETURN  60031 

end  oCC32 
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SUBRUUTINE  L)EL(rt.^)  OOCOl 

COMV.G.'t/SLVE  /  I  V  (  3C  )  .  TQL  (  2  .  )  .(,R  (  Od  I  2  )  ♦  3  (  “dS)  1  20  I  .X  (  12  1  !  tC  (  12  1  )  t  GO 0  02 

1  KL  (  26  )  tNAr;(  25*2)»A(0  6.121)»(MLt6';).vH65)tv2i63)*v3ib5)*  Oulu3 

2  V4(  65  )  iXMF  I  12  ^ )  *X  1  (  121  )  tX2  ;  121  1  ♦X-’  (  121  >  ♦Xl-AR  (  26  )  60.60  )  000  04 

Ei^-U  I  valence  (  I  V  (  1  :  .M  )  .  (  I  V  (  2  )  .f'  ENO  )  .  «  I  V  (  3  »  .iNCCP  )  •  (  I  V  I  4  )  .,\  .M  0  I  )  •  00'''00 

1  (lV(0)»NlT).(IV(6).NOTl.(lV(7).PF).(IV(8)»lli-K).lIVt9).Ii;'AX), 

2  (IVd  )  .  I  lRR'JR  »  ,  (  1  V{  1 1  )  .LASTCP  I  .  1  I'.M  12  )  .<:- )  DC007 

Dir-ENOION  /.(  1  )  .(j(  1)  DUC08 

DO  2L  J  =  l.u  D0CC9 

u*u01C 

DC  lu  I  =  l.M  D.Cll 

1  F  (  A  (  I  .  J  )  .NE  .0  .  )  /.‘.y  =  ....  +  A{  1  ,J  )  *  ^(  I  )  DC012 

lo  CGMINJE  DCC13 

W  (  J  )  =  '.V  .V  D  0  0  1  4 

2u  CONTI NJE  DOOlO 

RETURN  UOOlb 

END  00017 
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subroutine  RCALC  RJCOl 

COMKON/SLVE  /  I  V(  30)»TUL(2Q)tNR<5>5t2)t3(b5)»<N{12u)»X{121)tC(121)»  R0002 

1  K.LI  26)  .NAM(  2S.2)tA(S5,121).PIE(65).Vl(65).V2(6S)iV3(65).  R00C3 

2  V4(6S).XMF(120).Xl(121).X2n21)»X3(121)  tXuAR(  23  )  .R(  63  f63  )  RCC04 

EQUIVALENCE  (  I V ( 1 ) » M ) , { I V ( 2 ) »MENO )  .  ( I V ( 3 )  *NCCMP )t(IV(4)tN.NIUI)t  RCC05 

1  (IV(3).NIT)t(IV(6)*NOT).(IV(7)tPF).(IV(8)tITLR).(IV{9).I  Ti-iAX  )  »  ROC  06 

2  (  IV(  1C  )  ,  IlRROR  )  .  (  IV(  11  )  .LASTCP)  .  nv(  12  )  *<E  )  RUC07 

COMPUTE  R  R0C06 

DO  1  I  »  ItMEND  RCC09 

DO  2  J  «1.1  ROCIO 

RlItJ)“v./«0  RCOll 

2  continue  R0012 

1  CONTINUE  RC013 

DO  1C  K  « •  iNTOT  R0C14 

DO  11  1*1.K  R0013 

IF  (A(  I  .f.)  .EQ.O.  )  GO  TO  11  R0C16 

AUX  «  A(  I  ,K)  *  X(  <)  R0C17 

DO  12  J  *1.1  R0016 

IF  (  A  (  J  .K  )  .NE  .0.  )  Rll.Ji  =  A(J,IC)  *  AIF.X  +  R(l.J)  R0019 

12  CONTINUE  R0020 

11  continue  R0021 

1C  CONTINUE  RUC22 

DO  20  <  «  l.NCON'.P  R00  2  3 

IH  »  <  ♦  M  R0024 

MTA  *KL(<)  R0025 

MTB  »<L(<*1)  -  1  RQ026 

DO  21  L  «MTAtMTL  R0027 

DO  22  J  *1.M  R0028 

IF  ( A( J »L) .NE.O.  )  R(Ih.J)  =  R(Ih.J)  +  A(J»L)  *  X(L)  R0C29 

22  CONTINUE  RCC30 

21  CONTINUE  PC031 

20  CONTINUE  R0C32 

DO  30  J  »  2. MENU  R0033 

JL  =  J-1  R0034 

DO  31  I  »  l.JL  R0035 

R( I . J)  *  R( J. I )  R0036 

31  CONTINUE  R0037 

3w  CONTINUE  R0C36 

5c  RETURN  R0039 

END  R0040 
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SUBRQNTINE  CL0G(  a/.^ljAR  (  COCQl 

COMMON/ SLVE/  I  V(30)tTGL(20).NK(55.2).?(5‘)).KN(12U)»X(12n.C(121)*  CCC02 

1  KL  (  26  )  .NAM(  25*2)iA(55.121)tPlE(6!))*Vl«66)tV2l63).V3l6'j).  C0003 

2  VA(63).XMF(12w).Xl(121).X2(121).X3(121)«XoAR(23),R(65.65)  COCOA 

equivalence  (  I  V  (  1  )  tM)  I  (  I V  (  2  )  •MENU  )  .  (  I  V<  3  )  .NCO'-’P  ).(IV(A).N.NTuT),  CCC05 

1  (  IV(3)tNIT).(lV(6)»NOT)t(  IV(7).Pr).(lV(e>»ITLk).(  IVIVI.I TMAX )  t  C0C06 

2  (  I  V  (  iC  )  I  lEFRCR  )  .  (  1  V(  11  )  ♦LASTCP  )  » (  1  V(  12)  .i<;E)  C0007 

DlMLNolON  iV  (  1  )  ,i^’cAK(  1  )  ,AL?nA  (  1  )  CG008 

equivalence  .X2.ALPnA)  C0009 

DO  1  N  =  1  t  i'.CO.MP  CCCIO 

KLA  -  <L(N)  COCll 

KLB  =  <L(<+1)-1  C0012 

DO  2  J  *  KLA.^Lo  C0013 

ALPhAIJ)  =  C(J)  CuClA 

XXX  *  W(J)/aBAR(K)  ccai5 

I F ( XXX ,UT .C.C  )  ALPHA(J)  =  C ( J 1 >ALCG( XXX )  CCC16 

2  continue  CP017 

1  continue  CC018 

RETURN  C>.C19 

END  C-V.20 
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SUbKJUTINE  Lr^  C'C.NJ)  Luuul 

CO.V.Vo'.'./SLV:::/  lV(3C)tT<:L(2:it  .\a  (  0  j  .  2  )  .  ;  (  *;:)  12  .  )  .  X  (  12  1  )  *C  (  12  1  )  *  LOCO  2 

1  KL<  26  )  .NAM  (  25  *2  ).A(55.12n.PlLi6  3).Vl(6>).V2(o5).V3lt.5).  LOGO  3 

2  V4(65).XN;(-(12w).a1I121).X2(121).X3I121)  .x^AR(  2  3  )  .  ri  (  5  5 .6  3  )  LOO  04 

IMTtGLR  PF  L0n05 

ECU  1  valence  (  TUL  {  3  )  'X  -.in  )  .  (  TuL<  4  )  .XcTAPT  )  .  (  T'.'L  (  5  )  .oar  .;1N  )  L0006 

equivalence  (  IV(  1  )  (  IV(  2  1  .MENU)  .  (  I  V(  3  )  .NCCI'P  )  .  (  I  V  (  4  )  .N  ,■  NTOT  )  ,  LCC07 

1  (  1  V>  5  )  .NI  T  )  .  (  I  V  (  6  )  iNOT  )  .  (  1  V(  7  )  .PF  )  .  (  I  V(  0  )  •  ITlR  )  .  (  IV(9  )  .  I  TI-.AX  )  .  LCCOb 

2  ( IV( 1^5 . ILRRON) . ( IV( 11 )  .LASTCP)  .< !V( 12) .<E)  L0CC9 

DIMENSION  XX  (  1  )  .n>uJT  (  7  )  .CC  (  1 )  .P(  1  )  LUUlO 

ECU  I  valence (CC .XMF ).(Xa.X2).(P.V1)  LOO  11 

MON*  ^  L0012 

IF  (XSTART.LE.O.C)  XSTART  *  l.£-6  L0013 

DO  IG  I  »  l.M  L0C14 

P(  I  )  «  b(  I  )  L0C15 

A(  I  .NTOT4.1  )  =  C.O  LCri6 

DO  15  J  *  I.NTOT  L0G17 

A(I.NTOT+l)  X  A(I.NTOT  +  l)  +  Ad.J)  LCC18 

15  CONTINUE  LCC19 

lo  CONTINUE  L0C20 

DO  1  J  »  I.NTOT  L0C21 

CC(J)  *  C.^  L0022 

1  CONTINUE  L0C23 

CC(N+1)  «  -1.0  L0024 

C  ZEKO-Tm  simplex  is  to  uETERr.lNE  FEASIBILITY  L0023 

CALL  SIMPLE! -.M.N+1 . A . P . CC . NOU T . X X . P 1 E . V2 ♦ V 3 . V4 ♦ X 3 . R )  LOO 2 6 

ZT  *  XX(N-d)  L0027 

IFIPF.GE.C  IWRITE  (  NOT  .  1  06  )  XCU  T  (  2  )  .  Z  T  ,<OL'T  (  1  )  LOOZa 

lu6  FORMAT ( 12MCSIMPLEX  C.I4.29n  ITERATIONS.  MAX  MIN  ELEMENT* 1  PE  1 5 . B .  L0029 

1  12H.  CONDITION  .13)  L0030 

ZZT  *AMIN1 (ZT/2.0*  XSTART)  L0331 

DO  1C4  I  »  l.M  LOC32 

P(l)  *  P(I)  -  ZZT*A(1.N+1)  L0033 

104  CONTINUE  L0C34 

2'jj00201J*1.NTOT  L0035 

X(J)  .  XX(J)  LC036 

XMF(J)  X  l.B  L0C37 

2^1  CONTINUE  LCC36 

IF  (  Z  T  »  L  E  »  0  •  .  OR  .  K  UU  T  (  1  )  .  N£  *  </ )  00  TO  40  L0039 

C  SIMPLEX  LOOP  LG040 

FR2»1. E+20  L0C41 

DO  3b1  NN  *  1 .  NCuMP  L0042 

DO  3^2  J  *  1.  NTOT  LCG43 

CC(J)  *  C(J)  +  XMF(J)  -  l.C  L0C44 

3b2  continue  LCC45 

FN  =  FlOAT(NN)  -  1.0  L0C46 

CALL  SIMPLEII.M.N  .A.P.CC.NOuT.XX.PIE.VZ.VS.VA.Xj.R)  LG047 

IF  ( KOJT (  1  )  .NE.G  )  00  TO  50  L0046 

3uv  DO  3-3  J  =  l.NTUT  LU&49 

X(J)  »  XX(J)  LU050 

X(J)  =  (  FN*X1(J)  ♦  X(J)  )  /  (FN  l.D  L0051 

Xl(J)  *  X{J)  LGC52 

3u3  CONTINUE  LCC53 

call  RAR(X.XBAR)  LCC54 

K  =  1  LU055 

FR  -  ^.C  L0056 

DO  31G  J  =  l.N  .0057 

IF  ( J.OE.KLIX-*-!  )  )  <  =  X  1  .  -  L0058 

IF  (J.E«..',L(X).  ^Nu.  XoAiQ{  X).UT.C.r)FP  =  FK-XBAR(X)*.'LJG(Xi:AR(fx)  )  L0059 

IF  ( X ( J ) .OT . 0 .0 )  FR  =  FR  ♦  X FJ ) * ( ALGO ( X ( J )  )  ♦  C)J)  )  L0C60 
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XMF(J)  »  u.  LC061 

IF  (  XbAR(K)  .NE.O.  )  XMF(J)  =  X(J)  /  XbARIK)  L0062 

310  CONTINUE  L00o3 

IF  (PF.GE.O)  WRIT£(NOT.305)  NN . <0UT ( 2  )  . FR  L0C64 

FORMATIdH  S  I  VPLE  X  » I  3  » 1 H  t »  I A  1 12H  ITERATIONS  t8H  FR  ENC:=  1  PE  1 5 . 6  )  L006l> 

IF  (FR.GE.FR2)  GO  TO  399  L0066 

FR2=FR  L0067 

301  CONTINUE  L0068 

399  DO  ACC  J  »  1,N  L0069 

X(J)  X(J)  +  ZZT  L0070 

AUU  CONTINUE  L0G71 

RETURN  L0072 

U.  IF  (KOUT( 1 ) .GT.l )  GO  TO  50  L0073 

WRITE  (NOT.Al)  L0C7A 

^1  format ( 72H0THIS  PROBLEM  IS  INFEASIJLL.  THE  FOLLOWING  LINEAR  COMdI  L0075 
INATION  OF  RO.»S.  /IX)  L0076 

DO  lAu  I  xi,M  L0077 

IF  (PIE(  I  ) .NE.O. )  WRITE(NOT.lAl)  PI E ( I  )  .NR ( I . 1 )  .NR (  I .2  )  LOOTS 

lAl  FORMAT ( luX . 3H+  (.F15.8.5M  |  *  .2A6)  L0079 

lAJ  CONTINUE  L0080 

WRITE  (NOT.1A2)  LOOdl 

1A2  FuRMAT(AdHO  LEADS  TO  THE  FOLLOWING  INFEASitjLL  EQUATION.  /IX)  L00d2 

DO  15w  N  xl.NCOnP  LCC83 

MTA  =  KL(R)  LC08A 

MTa  »  XL(X-*-l)  -  1  L0085 

DO  151  J  »  MTA,  MTb  L0086 

D  *  0.  L0087 

DO  152  I  =1.M  00088 

D  =  PIt(I)»  A(I.J)  +  D  L0089 

152  CONTINUE  LC090 

IF  (O.NE.-.)  WRITE  (NOT.1A3)  0 , KN ( J ) ,NAM ( K . 1 ) . NAM ( K ,2 )  LCC91 

1A3  FORMAT { lOX . 3H+  (,F15.8.5N  )  *  .A6.AH  IN  ,2A6)  L0092 

151  continue  L0093 

15o  CONTINUE  L009A 

0*0.  LCC95 

DO  160  I  »1,M  L0C96 

D  «  PIE(  I  )»o(  1  )  D  L0097 

160  CONTINUE  L009d 

WRITE  (NOT.IAA)  D  L0099 

lAA  FORMAT! IHO, 15X,  7H+  0.0  ».F15.8)  LOlOO 

70  MON  *  1  LOlOl 

RETURN  L0102 

50  IF  (XOUT( 1 ) .NE.2)  GO  TO  60  L0103 

JT  *  XOUT ( 7 )  LOlOA 

DO  51  K  «  l.NCGMP  L0105 

IF  (  JT.GE.KLdC)  )  GO  TO  52  L0106 

51  CONTINUE  L0107 

52  WRITE  (NOT. 952)  KN ( U T ) , NAM ( X , 1 ) . NAM ( K  ,2  )  L0108 

952  FORMATdAH  THE  VARIAoLE  ,A6.AH  IN  »2A6,33H  IS  UNSOUNDED  AND  MUST  B  L0109 

IE  REMOVED)  LOllO 

GO  TO  7^  LOlll 

60  WRITE  (N0T,96v-)  L0112 

960  FORMAT(60H  SIMPLEX  ROUTINE  nAS  FAILED  DUE  TO  EXCESSIVE  ROUNO-OFF  E  L0113 


IRROR) 

GO  TO  70 

END 


LOllA 

L0115 

L0116 
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Calling  Sequence  for  Simplex  Subroutine 

The  simplex  subroutine,  SIMPLE,  may  be  used  to  solve 
a  general  linear  programming  problem  of  the  form:  Minimize 

n 

Z 

j=l 


subject  to 


n 

Z  ■  ‘’i  •  . m  (2) 

j-1 


The  a^j  is  stored  in  a  two-dimensional  array.  A,  with 
a.,  in  cell  A(i,j);  C.  is  stored  in  a  one-dimensional  array, 
C,  with  Cj  in  cell  C(j);  and  b^  is  stored  in  a  one¬ 
dimensional  array,  B,  with  b^  in  cell  B(i). 

The  calling  sequence  is 


CALL  SIMPLE(II,M,N,A,B,C,KO,X,P,JH,XX,Y,PE,E) 


where 


II  -  0; 

M  -  No.  of  rows,  m; 


N  -  No.  of  variables,  n; 
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A,  B,  C  Are  as  above; 

KO  =  A  subscripted  variable  of 
dimension  7; 

X  =  A  subscripted  variable  of  dimen¬ 
sion  n  or  more; 

P,  JH,  XX,  Y,  and  PE  =  Subscripted  variables  of 

dimension  m  or  more;  and 
E  =  A  subscripted  variable  of 
dimension  m  or  more. 

Upon  exiting  from  the  subroutine, 


X(l),X(2),...,X(n) 

P(l),P(2),...,P(m) 

K0(1) 


K0(2) 

K0(3) 


K0(4) 

K0(5) 


Contains  x^,X2,...,x^  (the  solution); 
Contains  the  shadow  prices; 

Contains  an  0  if  the  problem  was 
feasible,  1  if  the  problem  was 
infeasible,  2  if  the  problem  had 
an  infinite  solution,  and  3,  4,  or 
5  if  the  algorithm  did  not  terminate; 
The  number  of  iterations  taken; 

The  number  of  pivots  performed  since 
the  last  inversion; 

The  number  of  inversions  performed; 
The  number  of  pivot  steps  performed; 
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K0(6)  A  logicftl  variable  that  is  "true" 
if  and  only  if  the  problem  was 
feasible;  and 

K0(7)  Contains,  if  the  problem  had  an 
infinite  solution,  the  number  of 
the  variable  that  was  infinite. 

The  dimension  of  A  (line  X0009)  must  agree  (at  least 
in  the  first  subscript.)  with  the  dimension  of  A  in  the 
calling  program.  The  other  dimensions  need  not  agree  with 
those  of  the  calling  program. 

If  an  initial  basis  is  available,  this  basis  may  be 
communicated  to  the  subroutine  by  letting 


II 

X(i) 


1  , 

0.0  if  variable  i  is  not  in  basis, 
(non-zero)  if  variable  i  is  in  basis. 


and  the  other  quantities  remain  as  above. 

This  subroutine  differs  from  other  linear  programming 
routines  in  several  respects.  If  the  restraints  (2)  are 
linearly  dependent,  the  problem  is  considered  to  be  in¬ 
feasible.  This  is  the  case  because  the  chemical  equilibrium 
problem  cannot  be  solved  if  the  restraints  are  dependent. 

In  addition,  this  subroutine  was  written  to  be  as  scale-free 
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as  possible;  this  was  accomplished  by  computing  tolerances 
internally  in  the  subroutine. 
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C  AUTOVATIC  SIVPLEX  REDUNDANT  EQUATIONS  CAUSE  I NF EAS I J I L I T Y  XOOOl 

SUBROUTINE  SIMPLE! I NFLAG , MX .NN . A , 0 . C i KOUT i <B . P i JH , X i Y #PE * E )  X0002 

DIMENSION  ai 1 ) iCC 1) .KOUT t 71  I JHI  1  I iX{ 1 ) .P( I ) iY( 1 1 .  X0003 

1  Kb( 1 ) lEC 1  I iPEI 1  I iKOI 7)  X0004 

equivalence  (K»<0( 1  I ) . { ITER»<0( 2 ) I  I  (  INVCiKOI 3) ) .  X0C05 

2  (NUMVR.X0(4)  )  .(NUMPVi<0{5n»(FEAS.<0(6»  I  •(JT.<0(7)  )  X0006 

equivalence  (XX.LLI  xooot 

C  THE  FOLLOWING  DIMENSION  SHOULD  BE  THE  SAME  HERE  AS  IT  IS  IN  CALLER.  X0008 

DIMENSION  A(55.121)  X0009 

logical  FEAStVER.NEQ.TRIG.Xw.A3SC  XOOlO 

C  XOOll 

C  MOVE  INPUTS  ...  ZERO  OUTPUTS  X0012 

DO  1341  I  -  1.7  X0013 

K0(1)  >  C  XU014 

1341  CONTINUE  X0015 

M  ■  MX  X0016 

N  ■  NN  X0017 

TEXP  ■  .5**16  X0C18 

NCUT  ■  4*M  ♦  10  X0019 

NVER  ■  M/2  ♦  5  X0020 

M2  ■  M*#2  X0021 

IF  ( INFLAG. NE.C 1  GO  TO  1400  X0022 

C*  'NEW  START  PHASE  ONE  WiTn  SINGLETON  BASIS  X0023 

DO  14-2  J  ■  l.N  X0024 

KB(JI  >  w  X002S 

KQ  >  .false.  X0026 

DO  14w3  I  ■  l.M  X0027 

IF  ( A( I  .J) .EQ.0.0)  GO  TO  1403  X0028 

IF  (<Q.OR.A( I.JI.LT.O.O)  GO  TO  1402  X0029 

XQ  >  .TRUE.  X0030 

14u3  continue  X0031 

KB(JI  ■  1  X0032 

14w2  CONTINUE  X0033 

1400  IF  (  INFLAG. UT.l  )  GO  TO  1320  X0C34 

DO  1401  I  -l.M  X0035 

JH  ( I  I  ■  -1  X0036 

14ul  CONTINUE  X0037 

C*  'VER'  CREATE  INVERSE  FROM  • XB '  AND  'JH*  X0038 

132w  VER  •  .TRUE.  XC039 

1121  INVC  ■  0  X0040 

1122  NUMVR  ■  NUMVR  ♦!  X0041 

DO  11^1  I  ■  ltM2  X0042 

E( I )  >  O.w  X0043 

llwl  CONTINUE  X0044 

MM«1  X0045 

DO  1113  I  ■  l.M  X0046 

E(MM)  ■  l.v  X0047 

PE(l)  >  0.0  XU048 

X(l)  ■  Q(I)  X0049 

IF  (JH(I)  .NE.O)  JH(I)  ■  -1  X0050 

KM  •  MM  ♦  M  ♦  1  X0051 

1113  CONTINUE  X0052 

C  form  inverse  X0053 

DO  llu2  JT  •  l.N  X0054 

IF  (XUt JTI .EQ.O)  GO  TO  1102  X0055 

GO  TO  600  X0056 

C  60c  call  JMY  X0057 

C  CHOOSE  PIVOT  X005B 

1114  TY  ■  0.0  X0C59 

DO  1104  I  •  l.M  XCC60 
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IF  ( JH(  I  )  .NE.-l )  GO  TO  1104  X0061 

IF  {  ABS(Y(  I  )  ).LE.TY)  GO  TO  1104  X0062 

IR  s  I  XJC63 

TY  «  ABS ( Y ( I ) )  XC064 

11J4  C0N'‘NUE  XCC6b 

KB(JT)  «  0  X0066 

C  TEST  PIVOT  XuC67 

IF  (TY.LE.TPIV)  GO  TO  1102  X006a 

C  PIVOT  X0v69 

JH( IR)  »  JT  XC070 

Kb(JT)  =  IR  Xu071 

GO  TO  9v.O  XCC72 

C  9^^  CALL  PIV  XC073 

1K2  CONTINUE  X0C74 

C  RESET  artificials  X007S 

DO  11-^9  I  =  l.M  XCC76 

IF  ( JH( I ) . EG.-l  )  JM(I)  =  0  X0077 

liu9  CONTINUE  XU076 

•12^^  VER  *  .FALSE.  X0C79 

C  PcRrOKy.  ONE  ITERATION  XOOtiO 

C*  'XCK*  determine  FEASIolLlTY  X0081 

FEAS«  .TRUE.  X0&o2 

NEG  *  .FALSE.  X00ti3 

DO  12-1  I  =  liM  X0084 

IF  (X( I ) .LT.C.C)  go  TO  1250  X00d5 

IF  (JH(I).EQ.r)  FEAS  =  .false.  X0086 

12C1  CONTINUE  X0087 

C*  'GET*  GET  applicable  PRICES  X0C88 

IF  (.NOT. FEAS)  GO  TO  501  X0C89 

C  PRIMAL  PRICES  X0C90 

DO  503  I  *  l.M  X0091 

P( I )  «  PE(  I  )  X0092 

503  CONTINUE  X0093 

ABSC  =  .FALSE.  X0094 

GO  TO  599  X0095 

C  COMPOSITE  PRICES  X0096 

125u  FEAS  =  .false.  X0097 

NEG  =  .TRUE.  X0098 

501  DO  504  J  =  1,  M  X0099 

P(J)  =  0.  XOlOO 

5-4  CONTINUE  XOlOl 

AbSC  =  .TRUE.  X0102 

DO  5C5  1  =  l.M  X0103 

•'M  =  I  X0104 

IF  (X( I ) .OE.0.0)  GO  TO  507  X0105 

ABSC  a  .FALSE.  X0106 

DO  508  J  =  l.M  X0107 

P(J)  a  P(J)  +  E(MM)  XG108 

N'M  a  MM  +  M  X0109 

508  CONTINUE  XOllO 

GO  TO  5C5  XOl  1 1 

507  IF  ( JH( I ) .NE.O)  GO  TO  5C5  X0112 

IF  (X(I).NL.C.)  AoSC  a  .false.  X0113 

DO  510  J  =  1 .M  XOl  14 

P(J)  a  P(J)  -  E(MM)  XU115 

MM  a  M,M  ♦  M  X,j116 

51-  CONTINUE  X0117 

5-5  CONTINUE  XOllb 

C*  'MIN'  find  minimum  REDUCED  COST  X0119 

599  JT  a  0  X0120 
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dO  =  0.0  X0121 

DO  7C1  J  =1.N  X0122 

C  Skip  columns  in  basis  X0123 

IF  (Xb(J).Nt.C)  GO  TO  7C1  X0124 

OT  =  u.v.  X0125 

DO  3U3  I  =  l.V  X0126 

IF  ( A  (  I  t J  )  .NE.0.0  )  DT  «  DT  ♦  P(l)  *  A(I.J)  X0127 

3w3  continue  X0128 

IF  (FEAS)  uT  =  OT  C(J)  X0129 

IF  (ABSC)  DT  =  -  A65(DT)  X0130 

IF  (DT.GE.uo)  uO  TO  701  X0131 

bb=OT  X0132 

JT  =  J  X0133 

7ul  CONTINUE  X0134 

C  TEST  FOR  NO  PIVOT  COLUMN  X013S 

IF  (uT.LE.^)  go  to  203  X0136 

C  test  for  iteration  limit  exceeded  X0137 

IF  ( I TER.GE.NCUT )  GO  TO  160  X0138 

ITER  =  ITER  ■*■1  X0139 

C*  'UMY'  multiply  inverse  TImLs  A(..uT)  X0140 

6ww  DC  61^  I«  l.M  X0141 

Y ( I )  =  G.o  X0142 

61-  CONTINUE  X0143 

LL  =  V  X0144 

COST  =  C(UT)  X014S 

00  605  1=  l.M  X0146 

A1 JT  *  A(  I . JT  )  X0147 

IF  (AIJT.EG.O.)  GO  TO  602  X0148 

COST  s  COST  ♦  AIJT  *  PE(I)  XC149 

DO  606  J  »  l.M  X0150 

LL  =  LL  +  1  X0151 

Y(J)  s  Y(J»  ♦  AUT  •  l(LL)  X0152 

6^6  continue  X0153 

GO  TO  605  X0154 

602  LL  »  LL  ♦  M  X0155 

6u5  CONTINUE  X0156 

C  COMPUTE  PIVOT  TOLERANCE  X0157 

YMAX  =  0.0  X0158 

0*^  62'-  I  =  1  .fl  X0159 

YMAX  =  AMAXII  AuS ( Y (  I  )  I  .  YMAX  I  X0160 

62u  CONTINUE  X0161 

TPIV  =  YMAX  *  TEXP  X0162 

C  return  to  inversion  ROUTINE.  IF  INVERTING  XC.163 

IF  (VLR)  GO  TO  1114  X0164 

C  COST  tolerance  CONTROL  X0165 

IF  ( TRIG. AND. ub.GE. -TPIV)  GO  TO  203  X0166 

TRIG  =  .FALSE.  X0167 

IF  ( Bb.GE.-TP  I V )  TRIG  «  .TRUE.  X0168 

C*  'ROW  select  pivot  RCR  X0169 

C  AMuNG  EuS.  with  X=C.  find  maximum  Y  among  artificials,  or.  IF  NONE.  X0170 

C  GET  MAX  POSITIVE  Y(I)  AMONG  REALS.  X0171 

1000  IR  =  u  X0172 

AA  =  u.j  X0173 

KG  =  .FALSE.  X0174 

DO  1-50  I  =1.M  X0175 

IF  ( X (  I  ) .NE  .0 .O.OR. Y (  I  )  .LE.  TPI V  )  GO  TO  1050  X01  76 

IF  (UH(I).EO.O)  GO  TO  1044  X0177 

IF  (KQ)  GO  TO  1-50  X0178 

1245  IF  (Y(I).LL.AA)  GO  TO  1050  X0179 

GO  TO  lu47  X0180 
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1044  IF  (KC3)  GO  TO  1045  X0181 

KQ  *  .TRUE.  X0182 

1047  AA  =  Y( I)  X0183 

IR  =  I  X0184 

1050  CONTINUE  X0185 

IF  (IR.NE.O)  GO  TO  1099  X0186 

1001  AA  »  l.OE+20  X0187 

C  FIND  MIN.  PIVOT  AMONG  POSITIVE  EQUATIONS  X0188 

DO  1010  1  =  l.M  X0189 

IF  ( Y(  I  )  .LL.TPIV.OR.XI I ) .LE.O.O.OR.YI n*AA.LE.X( I )  )  GO  TO  1010  X0190 

AA  =  X ( 1  )  /Y { I )  X0191 

IR  »  I  X0192 

lOlo  CONTINUE  X0193 

IF  (.NOT.NEG)  GO  TO  1099  X0194 

C  FIND  PIVOT  AMONG  NEGATIVE  EQUATIONS.  IN  WHICH  X/Y  IS  LESS  THAN  THE  X0195 

C  MINIMUM  X/Y  IN  THE  POSITIVE  EQUATIONS.  THAT  HAS  THE  LARGEST  ABSF(Y)  X0196 

1016  Be  =  -  TPIV  X0197 

DO  1C3U  1  =  l.M  X0198 

IF  (X( I  )  .OE.O..OR.Y(  I) .uL.do.OR.YI n*AA.GT.X(  I )  )  GO  TO  1030  X0199 

Bb  =  Yin  X0200 

IR  »  I  X0201 

1030  CONTINUE  X0202 

C  TEST  FOR  NO  PIVOT  ROW  X020i' 

1099  IF  (IR.LE.o)  go  TO  207  X0204 

C*  ‘PIV*  PIVOT  ON  (IR.JT)  X0205 

C  leave  transformed  COLUMN  IN  Y(I)  X0206 

9C0  NUMPV  »  NUMPV  ♦  1  X0207 

YI  =  -Y( IR)  X0208 

Y( IR)  =  -l.u  X0209 

LL  «  0  X0210 

C  transform  inverse  X0211 

DO  9-4  J  =  l.M  X0212 

L  =  LI.  ♦  IR  X0213 

IF  (E(L) .NE.C.C)  GO  TO  905  X0214 

LL  =  LL  ♦  M  X0215 

GO  TO  904  X0216 

905  XY  =  E(L)  /  YI  X0217 

PE(J)  «  PE(J)  +  COST  *  XY  X0218 

E(L)  =  O.u  X0219 

DO  906  I  =  1  .M  X0220 

LL  »  LL  +  1  X0221 

E(LL)  =  EILL)  ♦  XY  •  Y(I)  X0222 

906  CONTINUE  X0223 

904  CONTINUE  X0224 

C  transform  X  X0225 

XY  =  X ( I R )  /  Y I  X0226 

DO  908  I  =  1,  M  X0227 

XNEW  =  X(I)  ♦  XY  *  Y(I)  X0228 

IF  ( VER.OR.XNEW.Ol. ...OR.Y(  I  ) .GT.TPIV.OR.XI  I )  .LT.O.  )  GO  TO  907  X0229 

X(  I  )  «  C.«>  X0230 

GO  TO  908  X0231 

9^7  X ( I )  «  XNEW  X02  32 

908  CONTINUE  X0233 

C  RESTORE  Y(IR)  X0234 

Y( IR)  =  -Y I  XU235 

X( IR)  *  -XY  X0236 

IF  (VER)  GO  TO  1102  X0237 

221  lA  «  JH( IR  )  X0238 

IF  (lA.GT.G)  KB(IA)  =  0  X0239 

213  KB(JT)  «  IR  X0240 
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JH( IR )  =  JT 

X0241 

IF  (  NU'-'PV.LE  )  GO  TO  1200 

X0242 

C  TEST  FCK  IWERSI^jN  ON  TriIS  ITCRATlv^N 

X0243 

INVC  =  INVC 

X0244 

IF  (  iNVC.E'i^.NVER)  GO  TO  1320 

X0245 

Tjvj  To  12^^^ 

X0246 

C*  d.NJ  Jf-  algor  I  Tn,v,  StT  EXIT  VAlgLo 

X0247 

C 

INFINITE  SOLUTION 

X0248 

2v7 

<  =  2 

X0249 

GO  TO  250 

X0250 

C 

PRJdLEM  13  cycling 

X025  1 

1 

K  =  4 

X0252 

GC  TO  25- 

X0253 

C 

FEASIBLE  -R  INFLASIJLl  SCLUTIuN 

X02  54 

2-3 

<  =  - 

X02  55 

25- 

IF  (.NOT.FEAS)  <  =  K  >  1 

X0256 

DO  1399  J  =  l.N 

X0257 

XX  = 

X0258 

<c)  J  =  Xl)  1  J  1 

X0259 

I  F  ( XBJ . Nl  .  -  1  XX  =  X ( XuJ ) 

X0260 

Xu ( J  )  =  LL 

X0261 

1399 

CONTINUE 

X0262 

C 

SET  'XOUT' 

X0263 

1  392 

DO  1393  I  =  1.7 

X0264 

XCUT ( I )  a  Xo( I ) 

X0265 

1393 

CONTINUE 

X0266 

RETURN 

X0267 

END 

X0?6fl 
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C  yATKIX  IiMVEKSlON  aITi-i  ACC  jM  f"ANY  1  No  f.CLCTIuN  OF  LINEAR  EoJATlONS  MOOOl 

SUdROUT  INE  MAT  INV  (  A  f  N  .  b  . X.  •  I  NA  .  I Nd  .  I  P  .  I  5  1  NO  )  X.''yr02 

C  I*.  C  C  0  3 

DIMENSION  B( 1 ) t INA( 1  I . IND( 1  I . IP(  1  I  XOCDA 

logical  ip  XOCOS 

DIMENSION  A(6&,65)  .•luCOS 

C  MOO  3 7 

C  initialization  MOOOb 

DO  20  J  *  1»N  W0009 

iP(j)  =  .false.  moo  10 

2u  CONTINUE  XOOll 

C  bIG  LOOP  UN  1  Mu:i2 

DO  bib  I  =  l.N  M0013 

AMAX  X  J.„  MOOIA 

C  SEARCn  FOR  PIVOT  ELEMENT  MOOIS 

DO  lu5  J  =  liN  MC 016 

IF  ( IP( J)  I  oO  TO  1-3  M0C17 

DO  1-0  X  =  l.N  MOCia 

IF  (IP(K)  .OR.  AbS  ( AMAX  )  .  GE./E':  ( A(  J  .K)  )  )  GO  TO  ICO  M0019 

IROvy  =  J  M0020 

I  COL  =  X  MOO  21 

AMAX  =  A(J.XI  M0C22 

lOo  CONTINUE  M0023 

1j5  continue  M0024 

IF  (AMAX. c(w. 0.0)  GO  TO  76c  M0025 

IP( ICOL)  =  .TRUE.  M0026 

C  INTFRChANGE  rows  to  PUT  PIVOT  ELEMENT  CN  DIAGONAL  M0027 

r'  (  IROW.EU.ICOL)  GO  TO  26G  MOOZfl 

'lO  200  L  =  l.N  M0029 

SWAP  X  A( IROW.L )  M0030 

A( IROW.L)  X  A( ICOL.L)  M0031 

A (  ICOL  tL  )  X  SWAP  M0032 

200  CONTINUE  M0033 

IF  (M.EQ.O)  GO  TO  260  M003A 

SWAP  X  8  (  IR'Jrt  )  M0035 

b( IROW)  *  B( ICOL)  M0036 

b(ICOL)  =  SWAP  M0037 

26o  INA(I)  X  IRUW  M0038 

INb(I)  X  ICOL  M0037 

C  DIVIDE  PIVOT  ROW  bY  PIVOT  ELEMENT  M0040 

A (  ICOL .ICOL)  X  1 .0  M0041 

DO  36-  L  =  l.N  MCC42 

A(  ICOL.L)  X  A(ICCL.L)  /  AMAX  MCC43 

30-  CONTINUE  MOO 44 

IF  (  I.NE.-)  b(ICUL)  =  3(ICCL)  /  A*'AX  M0C46 

C  COXPL'ITE  Tm^  PIVOT  M0C46 

3d0  DC  660  LL  =  l.N  M0C47 

IF  (LL.^-.IOCL)  oU  6  60  MC04o 

£.■.4?  X  A(LL.ICOL)  M0045 

A(LL.ICUL)  =  .ijCSO 

--43^  L  =  1.N  .‘'■.0061 

A(LL.L)  X  A(LL.L)  -  A(IC-L.L)  "  S..AP  M0062 

46-  CoNTINJl  M0063 

IF  (M.NE.-)  b(LL)  =  (  LL  )  -  3  (ICOL)  *  S.-.A'’  M0C64 

66-  CONTINUE  “.-066 

0  76  CONTINUE  .••■.CC66 

6--  IF  C'.LT.-)  RET'jRs  MOO  6  7 

C  interchange  C-LU'-'NS  M006b 

DC  71-  I  X  l.N  .M-'l69 
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L  =  N  ♦  1  -  1 

IF  (  INA(L)  .tQ.INB(L)  )  OU  TO  710 

IRC-v  =  I\A(L)  ■  N0u62 

ICCL  =  INb(L)  MCC63 

DO  705  <  =  l.N  y.0C64 

S^AP  =  A(k,,IROw)  MCC65 

A  (  K  ,  1  RC»^-)  =  A(<tICOL)  M0G66 

A(K,1C0L)  =  SWAP  M0067 

7.5  continue  y,006t} 

71^  CONTINUE  M0069 

74 j  RETURN  y007C 

c  singularity  flag  M0071 

75w  ISING  =  1  ♦N  -  I  i'lu0  72 

G  C  T  0  6  .  G  y  0  0  7  3 

END  M0074 
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Appendix  B 

MATRIX  NOTATION  AND  FURTHER  PROOFS 


The  derivations  in  the  preceding  sections  would  be 
facilitated  by  the  use  of  matrix  notation  rather  than  sub¬ 
scripted  variables.  We  introduce  the  following  symbols  to 
correspond  to  the  subscripted  variables  used  in  Sec.  3. 


Subscripted  Variable 

Matrix 

Size  of  Matrix 

"ij 

A 

MxN 

•>1 

B 

Mxl 

’'j 

Y 

Nxl 

D 

Nxl 

C 

Nxl 

ir 

Mxl 

R 

MxM 

X. 

J 

X 

Nxl 

The  single-column  matrices  may  also  be  thought  of  as  vectors. 
We  use  here  the  convention  that  an  operator  applied  to  a 
matrix  means  that  the  operator  operates  on  each  element  of 
the  matrix.  For  example,  log  Y  is  the  Nxl  matrix  consist¬ 
ing  of 
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/  \ 
log  Yi 

log  y2 

log 

\ 

T 

The  superscript  indicates  the  transposition  of  a  matrix. 
We  assume  that  the  elementary  results  of  matrix  theory  are 
known.  For  example,  it  is  known  that  the  inverse  of  an 
invertable  symmetric  matrix  is  symmetric.  The  square 
diagonal  matrix  whose  diagonal  is  one  of  the  vectors  pre¬ 
viously  defined  will  be  denoted  by  the  previously  defined 
vector  in  elongated  type;  that  is, 

D  =  dlag  (D) 

and 

Y  =  dlag  (Y) 

Equations  (3.2)  and  (3.7)  in  matrix  notation  are 


AX  -  B 


(B.l) 


X  -  Y  (D'^A^’n  -O'^c  -  log  Y)  . 


(B.2) 


To  see  the  ease  of  matrix  notation,  we  may  substitute  (B.2) 
into  (B. 1)  to  get 


AY0"^A^ff  -  B  +  aY(0'^C  +  log  Y)  . 


By  letting 


R  -  aYO'^a^ 


and 


S  -  B  +  aY(D'^C  +  log  Y)  , 


we  see  that 


Rtt  =  S 


(B.3) 


(B.4) 


(B.5) 


(B.6) 


corresponds  to  (3.10). 

In  Sec.  4,  we  evaluated 
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but  we  did  not  give  the  details  of  the  computation.  The 
algebra  of  this  evaluation  is  very  difficult  unless  matrix 
algebra  is  used.  In  matrix  notation,  (B.7)  is  O^DY  ^0, 
where  0  =  X-Y.  From  (B.2)  we  have 


0  =  Y  (D'V"  -D'^c  -  log  Y)  -  Y  . 


(B.8) 


Hence , 

0’'OY'^e  =  (/a0'‘  -  c^O’^  -  log  y’’)YDY'^9  -  y^OY'^q 

=  iT’’A(D'^YDY'be  -  (c’'d'^  +  log  Y’’)DYY'^e  -  y^’y^’-Qs 

=  it’’ AO  -  +  log  y’’)D0  -  □’’o  .  (B.9) 

Since  AX  =  B,  A0  =  AX-AY  =  B-AY.  Also,  in  the  chemical 
equilibrium  formulation, 

j=l  j=n+l  k=l'je(k^  / 


and 
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Hence, 


(c’'d''-  +  log  y’')D9 


n 


N 


=  ^  (Cj  +  log  y^)0^  +  ^  log  yj(-Qj) 

j=l  j=n+l 


k=l  \jc^k) 


+  log  y^)  -  0^  log 


Z  (  Z  ®k> 

k=l  \jf(k) 


n 


z 

j=l 


0 . (c .  +  log  y . ) 
J  J  J 


m 


n 


n 


=  z 


TT  .  (b  . 
11  1 


Z^/j)-  z 


0^(c^  +  log  y^)  (B.IO) 


j  =  l  i=l  \  j=l 


j  =  l 


in  the  context  of  the  chemical  equilibrium  problem  used  in 
Sec.  4. 


Next  we  wish  to  show  that 
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“  e?d. 

0 

L  y. 
j=l  J 

as  stated  in  (4.14).  First,  we 

prove 

Lemma  1:  Let  y^,yo 

.  •  •  • ,  y^- 

be  positive  numbers  and  let 

0^,02,...,0^  be  any  real 

n  umbers 

.  Let 

^  0^  ( 

r  -  V  — 1.  -i 

y  0.1 
i-i  V 

"  -  2,  yj  - 

j=i  J 

r 

j=l  J 

• 

Then , 

1)  G  ?•  0 

ii)  G  =  0  if  and 

only  if 

_  «2 

0 

r 

yi  yj 

•  •  • 

"  ^r  ' 

Proof:  Let  a.  =  0./y., 
-  J  J  J 

j=l,2,.. 

. ,r .  Then, 

/  "" 

\2 

r 

r  -  V 

I  7  Of .y . 

Vl=i  ' 

)_ 

L  2  3 
j=l 

r 

j»l  ^ 
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-1 


I 


Li=l 


-1 


I 


Li=l 


lot.a.y.y.  +  a.y.y. 
1  i-'i-'j 


which  is  result  i) .  The  proof  is  completed  by  noting  that 
G  =  0  if  and  only  if  for  all  i  and  j;  this  proves 

ii)  . 

Now  we  can  prove 

Theorem  1;  In  the  chemical  equilibrium  problem 


i) 


y-i_i 

L  y. 
j=l  J 


0 


^  0^d. 


=  0 


j  =  l 


if  and  only  if  there  exist 


ii) 
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numbers  Of,  ,  Wo  , .  .  .  ,or  such  that 

L  z  p 


a)  0.  =  Wr .  -,y.  j^n 

J  Ij]  J 


b)  0 .  =  w ,  S .  .  j>n 

J  J-n  j-n 


Proof :  The  proof  follows  by  noting  that  for  i  >  n 


0. 

I 


I 


Then , 


by  lemma  1.  Furthermore,  by  lemma  1,  if  the  equality  holds, 

then  for  each  k  there  is  a  number  a.  such  that  0.  =  a,  y .  if 

k  J  k-'j 

j  c  k.  This,  noting  that  b)  follows  from  the  fact  that 
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0.  *  y  0.  for  i>n  , 

^  L  2 

jc<i-n) 

completes  the  proof  of  the  theorem. 

Our  final  result  is 

Theorem  2:  In  the  chemical  equilibrium  problem,  with 

(yT>yo>*-->y  )  feasible  and  0, ,0«,...,0  calculated  as  in 
•'I  -^2  '■'n  12  n 

(4.7) 

n 

1)  ^  (c^  +  log  y^)  i  0 

j'l 

n 

ii)  ^  ® j  j  ^  ^  only  if 

j-1 

(y^,y2>  •  •  •  ’Yr,)  optimal. 

Proof:  i)  follows  from  Theorem  1,  (B.IO),  and  the  fact 

that  (y^^ ,y2 >  •  •  •  »yj^)  fs  feasible. 

To  prove  ii) ,  we  assume  that 

n 

^  0^(c^  +  log  y^)  -  0  . 

j-1 


Then , 
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N 

I 

j=l 


0^d. 
J  J 


0  . 


and  is  as  in  ii)  of  Theorem  1.  Combining  b)  of  Theorem 
1  and  (4.12)  we  have 


0 


k+n 


^k^nri-k 


°  “k^k 


or 


"k 


v 


m+k 


Next,  we  combine  a)  of  Theorem  1  with  (4.7)  to  get 


y  a..  -  c.  -  log  y.  +  ttI,,. 

L  1  ij  J  ^  2  [j]+m 

i=l 


or 


I 


ttI  a  . 
1-  ij 


i-1 


c . 
J 


1  ^ 
log  y, 


0  . 
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Thls  last  result  is  the  optimality  condition  for 
‘ * ’^n^  as  given  by  (1.4),  and  this  demonstrates 
the  forward  implication  of  ii).  The  converse  follows  from 
the  fact  that  optimality  implies  that  the  objective  function 
cannot  be  decreased. 
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